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ABSTRACT 


A state-of-the-art computer code has been developed that incorporates a modified 
Runge-Kutta time integration scheme. Upwind numerical techniques. Multigrid acceler- 
ation, and Multi-block capabilities (RUMM). A three-dimensional thin-layer formulation 
of the Navier-Stokes equations is employed. For turbulent flow cases, the Baldwin-Lomax 
algebraic turbulence model is used. Two different upwind techniques are available, van 
Leer’s flux-vector splitting and Roe’s flux-difference splitting. Full approximation multi- 
grid plus implicit residual and corrector smoothing were implemented to enhance the rate 
of convergence. Multi-block capabilities were developed to provide geometric flexibil- 
ity. This feature allows the developed computer code to accommodate any grid topology 
or grid configuration with multiple topologies. The results shown in this dissertation 
were chosen to validate the computer code and display is geometric flexibility, which is 


provided by the multi-block structure. 
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Chapter 1 INTRODUCTION 

1.1 Historical Perspective 

Over the past thirty to forty years the advancement in computer resources and capa- 
bilities has grown exponentially. Along with this advancement has been the development 
and implementation of many numerical techniques for predicting fluid flows for different 
geometrical configurations. As both hardware and software improve, more and more 
complicated problems are being analyzed. Aerodynamic flows ranging from continuum 
to rarefied flows are being simulated with computers. One of the first major developments 
for continuum flow in computational aerodynamics was the boundary integral method, 
also known as the panel method. Its initial use was for the solution of subsonic lin- 
earized potential flows. This method was first employed by Hess and Smith in 1962 [1] 
for computing flows for three-dimensional, non-lifting bodies. Panel methods were ex- 
tended to lifting flows [2] for inviscid, low Mach numbers, where compressibility effects 
were small, and supersonic flows [3] (presently there exists integral equation methods 
for solving transonic flows, such as those by Kandil and Yates [4] and Kandil and Hong 
[5]). Transonic flows presented difficulties because the subsonic flow regions required 
elliptic solution techniques, and the supersonic flow regions required hyperbolic solution 
techniques. 

In 1970, Murman and Cole [6] successfully solved the transonic small disturbance 
equations for transonic aerodynamic flow fields using a successive line over relaxation 
(SLOR) algorithm, with a scheme known as the “type difference scheme . They used 
central differencing for the subsonic regions and upwind, one sided differencing, for the 
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supersonic regions. This provided the proper biasing for the numerical differentiation of 
the discretized governing equations relative to the characteristic directions for information 
propagation. For many flow configurations the potential flow equations are sufficient, but 
they assume isentropic, irrotational flow which is not valid for flows containing shock 
waves. However, for weak shock waves, these assumptions are valid up to second-order 
approximations. A more precise solution of inviscid transonic flows can be obtained by 
using the Euler equations. 

Solving the Euler equations requires more memory and computational time than 
the potential equations. Another landmark in 1970 was the work of Magnus and 
Yoshihara [7], who produced one of the first Euler computer codes accepted for computing 
transonic flows. They used the Lax-Wendroff scheme, which required an added artificial 
viscosity to remain stable. Artificial viscosity, also called artificial dissipation and 
numerical dissipation, is a numerical term that is related to the type of technique used 
in approximating the governing equations of motion, and should not be mistaken with 
the physical viscosity of a fluid. It will be referred to as numerical dissipation for the 
remainder of this dissertation. Numerical dissipation is required to stabilize numerical 
schemes, and is not a physical phenomena of the flow field. Another significant 
contribution came from MacCormack in 1969 [8], where he introduced a two stage 
predictor-corrector explicit scheme for iteratively solving inviscid flows about three- 
dimensional bodies. In 1972, MacCormack and Paullay [9] developed the rationale 
used for applying space discretization to the Euler equations. Another major contribution 
came from Beam and Warming [10], who used an implicit finite-difference algorithm to 
solve the conservative form of the Euler equations. A trapezoidal formula was chosen 
to integrate the unknown conserved variables, which produced an implicit difference 
equation. An Alternating Direction Implicit (ADI) method, based on those introduced by 
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Douglas [11], Peaceman and Rachford [12], and Douglas and Gunn [13], was used to 
execute the integration in time. Beam and Warming [10] incorporated a hybrid scheme 
which switched from central differencing, for subsonic regions, to upwind differencing 
whenever the local characteristic speeds were of the same sign, as is the case for 
supersonic regions. This is similar to the technique used by Murman and Cole [6]. 
The reduction of an implicit scheme into an alternating direction scheme was originally 
introduced by Gourlay and Mitchell in 1966 [14]. Briley and MacDonald [15] developed 
an equivalent alternating direction technique for solving non-linear hyperbolic equations, 
and applied it to the three-dimensional, compressible Navier-Stokes equations in 1974. 
They applied central differences to compute the spatial flux derivatives. The alternating 
direction method was also used by Steger [16] in 1978. He incorporated it into general 
curvilinear coordinates for computing transonic flow about arbitrary two-dimensional 
geometries using a finite difference scheme. It was also employed by Pulliam and Steger 
[17] in the same year for computing transonic, three-dimensional inviscid and viscous 
flows using a finite-difference method, with central differencing applied to the spatial flux 
derivatives. Beam and Warming [18] again used the alternating direction method in their 
finite difference scheme for computing the solution to the compressible Navier-Stokes 
equations, where central differencing was applied to the spatial flux derivatives. 

Expanding on the concept of biasing the type of numerical differencing used on 
the flux vectors, were the biasing depends upon whether the flow field surrounding the 
point of interest is subsonic or supersonic, Steger [19] introduced the concept of splitting 
each flux-vector, from the conservation law form of the governing equations, into two 
flux-vectors. The vectors were chosen so that the Jacobian matrix of one of the split 
flux-vectors would contain only positive real eigenvalues, and the other only negative 
real eigenvalues. This type of separation allows for upwind differencing to be used for 
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each point of interest: backward differences would be used for terms associated with 
positive eigenvalues and forward differences would be used for terms associated with 
negative eigenvalues. Executing this type of splitting of the flux-vectors removes the 
need to switch between central differencing and upwind differencing, where the type of 
differencing chosen was based upon the local Mach number at the point of interest. Steger 
and Warming [20] developed this technique in what is known as flux-vector splitting in 
their paper in 1981. Upwind differencing is an attempt to model the directions of signal 
propagation. Two schemes that closely model the characteristic propagation directions 
are the A-scheme, by Moretti [21], and the Split Coefficient Matrix (SCM) scheme, by 
Chakravarthy, Anderson, and Salas [22], Unfortunately, these two schemes can only be 
applied to the non-conservative form of the governing equations, and therefore require 
shock fitting techniques to locate the shock waves. In 1982, van Leer [23] introduced 
another type of flux-vector splitting that provided smooth transitions between the split 
fluxes when the eigenvalues changed signs, and good shock capturing capabilities. This 
approach is considered a pseudo particle approach. Other types of upwind methods 
are those that iteratively solve the Riemann problem, such as Godunov [24] proposed 
in 1959. This approach is based on the shock tube membrane rupturing problem, and 
solves the Riemann problem at every cell face of the physical domain, searching for 
a shock, expansion, and7or contact wave. Another approach is to approximate the 
Riemann problem to a set of equations that can be solved exactly. Schemes that solve the 
approximate Riemann problem are known as flux-difference splitting schemes. Examples 
of these types of methods have been developed by Roe [25, 26], Lombard, Oliger, and 
Yang [27], and Engquist and Osher [28], Both flux-vector splitting and flux-difference 
splitting can be applied to the conservative form of the governing equations, and therefore. 
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as shown by Lax [29], Lax and Richtmyer [30], and Lax and Wendroff [31], capture a 
shock implicitly by solving the governing equations. 

A method introduced in 1981 by Jameson, Schmidt, and Turkel [32], applies central 
differencing to the spatial flux derivatives and used the classical four-stage Runge-Kutta 
time-stepping scheme for solving the governing equations. This approach was modified to 
require less computer memory, as shown by Jameson and Baker [33], with the knowledge 
that the coefficients for each time stage could be modified to provide various stability 
and amplification characteristics. In 1985, Jameson [34] explained that adjusting the time 
integration coefficients and the number of stages would alter the stability and amplification 
regions. He also revealed the benefits of evaluating the numerical dissipation at various 
stages, using a different set of coefficients than that of the time integration. This work 
provided a significant step for the central difference computer programs. The use of 
multistage Runge-Kutta methods with modified time integration coefficients has also been 
investigated for upwind solvers [35-37], The flux-vector splitting and flux-difference 
splitting methods are becoming more popular and the multistage Runge-Kutta scheme 
is widely accepted. 

The solution of realistic fluid dynamic problems can become CPU intensive, and as 
more grid points or cells are added, the amount of CPU time increases in a non-linear 
fashion. One approach being used to accelerate the convergence rate of iterative schemes 
was introduced in 1964 by Fedorenko [38]. He presented a technique called multigrid. 
This process was further developed by Brandt [39] for boundary value problems and 
applied to the small disturbance equation by South and Brandt [40] for transonic flow 
calculations. It was later applied to the Euler equations by Ni [41] and Jameson [34, 42]. 
This method has become an integral part of many steady-state flow solvers for Euler and 
Navier-Stokes equations, for both central-difference and upwind-differencing schemes. 
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As the numerical techniques advanced beyond the geometrically uncomplicated test 
configurations to multiple element airfoils, intemal/extemal engine designs, and entire 
aircraft configurations, the grid generation process became much more difficult. Much 
research has been directed toward this problem [43-46], It became apparent that 
developing the grid topologies for these configurations would be easier if the geometries 
were divided into sections or blocks and then joined together in a fashion that a computer 
program could analyze. This approach is often referred to as domain decomposition. 

There are different types of domain decomposition. One approach is to have the grid 
patches or blocks overlap and/or be embedded with each other. This is also referred to as 
the Chimera grid scheme. One of the first to use this technique was Boppe in 1977 [47] 
for transonic wing flows. Others were Hedman [48] and Thompson[49]. Atta [50] and 
Atta and Vadyak [51] used this approach to solve the transonic full potential equations. 
Benek, Steger, and Dougherty [52] and Benek, Buning, and Steger [53] used the Chimera 
approach with the Euler equations for transonic airfoil and wing/body configurations, 
respectively. Eberhardt and Baganoff [54] used the embedded grid approach for a 
supersonic blunt cylindrical body configuration, and tried to address the problem of 
properly maintaining flow discontinuities across grid boundary interfaces. Maintaining 
conservation across grid interface boundaries is a difficult problem for the Chimera grid 
embedding scheme, especially for higher order extrapolations. 

Another type of domain decomposition is grid or block patching. This approach 
does not allow the grid patches or blocks to overlap. The blocks interface along the same 
surface. Grid or block patching was done by Chambier, Ghazzi, Veuillot, and Viviand 
[55] in 1981, for a system of hyperbolic equations. They used compatibility equations 
to develop the interface boundary approach, which provided good results for transonic 
channel flows. Unfortunately this approach is not conservative and therefore unsuitable 
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where flow discontinuities cross block interfaces. This problem was approached by 
Rai, Chakravarthy, and Hessenius [56] and Rai [57]. Mastin and McConnaughey [58] 
investigated the issue of higher order solutions on block patched grids with C 1 continuity 
(lines meet one to one) at the interface with no grid discontinuities and compared these 
results to configurations where the blocks overlapped. A number of researchers have 
used the block patching method without the C 1 continuity condition [59-66]. Others, 
while still maintaining the C 1 continuity were able to handle a variety of complex 
configurations [67-70], It should be noted that only C 1 continuity provides both higher- 
order extrapolations and conservation of fluxes across block interfaces. 


1.2 Physical Problems of Interest 

Many present day Computational Fluid Dynamics (CFD) computer codes have state- 
of-the-art solvers. The current effort being put forth is to take these numerically advanced 
computer codes and use them on more physically demanding geometries. It is no longer 
sufficient for a computer code to only provide analysis for a wing; it needs to be capable 
of including a fuselage and a nacelle. In analyzing such a configuration, one can see 
that these geometries cannot be accommodated with only one grid or mesh. In many 
cases, different grid topologies should be used for different parts of the configuration. 
For example, a C-0 mesh may be desired for the wing, an H-0 mesh for the fuselage 
and a polar grid inside the nacelle. To accommodate these different components, even 
if they were of the same mesh topology, a computer code capable of handling the 
different sections of the configuration separately with sufficient communication between 
the sections is required. This computer code must either be specially designed for this 
particular type of problem, or be what is generally called a multi-block computer code. 
A multi-block computer code would be the most accommodating. It should be flexible 
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enough to allow the user to divide a given configuration into many different sections, or 
blocks, and use a mesh topology that best suits each particular section or block. Plus, 
if required, each block could be handled differently, in terms of the flow solvers or 
governing equations. This approach is often referred to as domain decomposition. 

Another configuration that requires multiple block capabilities is afterbodies, with 
internal nozzles. These configurations also require different mesh topologies. For the 
external body a polar grid is used. For the internal nozzle, a polar grid is best suited at 
the wall of the body for two dependent reasons. One is that at the end of the body, where 
the external and internal blocks meet, the meshes need to match with C 1 continuity, due 
to the interface constraints of the current multi-block computer code. Therefore, the best 
way for the grid lines to match up at interfaces is to use the same mesh topology, and due 
to the external geometry constraints, a polar grid is best suited. Unfortunately a polar grid 
in the internal nozzle is extremely demanding because the nozzle is not axisymmetric, but 
rectangular; therefore two mesh topologies are used in the internal nozzle. A polar mesh 
is used at the wall region, because it will match the external geometry with C 1 continuity 
and it requires packing in only one direction. This is beneficial in using Navier-Stokes 
equations, especially if using an algebraic turbulence model, because there will be only 
one coordinate direction necessary for a length scale. The second mesh topology is an 
H-H mesh, which will interface with the polar grid and fill in the remainder ol the interior 
of the nozzle. This case will be explained in more detail in the results section. 

Other configurations of current interest, such as the National Aerospace Plane, with its 
multiple scram jets, and Advanced Tactical Fighters with thrust vectoring/thrust reversing 
nozzles definitely require a multi-block computer code for computational analysis. In 
analyzing scram jets and thrust vectoring/thrust reversing nozzles, a multi-block computer 
code is needed to handle the different grids and the multiple boundary conditions a 
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particular block face might encounter. It would not be uncommon for a block face to 
have two or more different boundary conditions. 

1.3 Objective of Dissertation 

The objective of this dissertation was to develop a state-of-the-art computer code 
that was capable of handling all of the aforementioned configurations. The efforts were 
directed toward steady-state solutions, which allowed for different types of convergence 
accelerators to be used. The main thrust of this project was to develop a computer 
code that was capable of handling many different problems of various mesh topologies, 
configurations, and boundary conditions, without requiring any changes to the basic 
computer code. To achieve this, the computer code had to have grid independent 
subroutines that were adaptable to various boundary conditions. It had to allow for more 
than one type of boundary condition on a given grid surface, such as for the surface ol 
the grid wrapping around a wing and forming a wake line or wake cut. Mesh topology 
independence pertains not only to one type of mesh at a time, but also the ability to 
handle multiple mesh topologies at the same time. To accommodate all of these desired 
qualities it was determined that a multi-block computer code was required. A multi-block 
computer code provides the flexibility of handling all of the different grid configurations, 
plus the interaction of cases that if not analyzed with such a computer code, would require 
a computer program modified just for the one specific configuration. Also, a multi-block 
computer code can provide the flexibility of having different mesh topologies interact. 
Knowing that this type of computer code is going to be used provides more flexibility for 
the grid generation process, by allowing the best grid topology for each particular section 
of the configuration being studied, without being overly concerned with what topologies 
the other grid sections are going to posses. This remains true even for a multi-block 


9 


computer code that requires C 1 continuity at the block interfaces. C 1 continuity does not 
put a restriction on the types of meshes that can be used. 

1.4 Requirements of Proposed Computational 
Fluid Dynamics Computer Code 

The proposed Computational Fluid Dynamics (CFD) computer code was to be capable 
of accommodating intemal/extemal flows, wing body configurations, and afterbodies with 
internal nozzles. This computer code was to be mesh topology independent, allowing it to 
handle C-O, C-H, O-O, H-H, and H-0 mesh topologies and their interactions. Although 
today’s state-of-the-art computers permit very large memory requirements, a judicious 
use of computer memory is still required. As more complex configurations are tested, it 
is easy to have a half million to one and a half million points in a grid. In developing 
this computer program consideration was given to its readability. 

There is a trade off at some point between legibility and computational efficiency. An 
example of this is in having a three-dimensional array in a corresponding set of nested do 
loops”. If the array is kept as a three-dimensional array, then when it is compiled only 
the inner “do loop” will be vectorized, but if the three-dimensional array is collapsed into 
a one-dimensional array then only one “do loop” is needed and instead of vectorizing 
a line at a time, the computer will vectorize a volume at a time. This would greatly 
decrease the amount of CPU time required per calculation. Unfortunately working with 
and manipulating a collapsed three-dimensional array can become quite cumbersome. The 
method adopted in the present work was to maintain the arrays in their three-dimensional 
form in the subroutines, attain vectorized inner “do loops” where ever possible, and write 
the computer code in such a manner that would allow the rearranging of the arrays in a 
subroutine from three-dimensional to either two-dimensional or one-dimensional, without 
affecting the main program or the other subroutines. 
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1.5 Numerical Analysis 

Much has been accomplished in using the central difference operators with explicitly 
added dissipation for the solution of transonic flows, but as with all approaches there 
are a few drawbacks. One is that central difference computer codes require added 
second and fourth order dissipation for stability and to reduce oscillations. This requires 
a certain amount of numerical experimentation for tuning the coefficients. Another 
common problem is that due to the numerical operator’s nature, it tends to smear contact 
discontinuities and shocks. The central-difference operators generally require more points 
at and in the shock so as to produce a sharp shock. A significant advantage of the 
central-difference methods is that they require less logic in evaluating the fluxes than 
the current upwind solvers; where upwind is defined as using only one-sided differences. 
Also, central-difference operators are generally set up with enough damping to provide 
sufficient smoothing, even in an explicit time integration approach. Thus, they again can 
save CPU effort by not requiring an implicit method just to damp the high frequency 
errors generated during the iteration process. 

Some of the advantages of an upwind scheme are that its damping characteristics are 
built into the flux evaluator. There is no need to add explicit second and fourth order 
dissipation for stability; therefore there is no fine tuning of the dissipation. Also, by 
evaluating the fluxes with an upwind approach there is less smearing of shocks and contact 
discontinuities, and fewer points, as few as two or three, are required to adequately capture 
a shock than are required for a central-difference scheme. This can be very helpful in 
making calculations on a preliminary grid which is constructed without prior knowledge 
of where certain physical changes of the flow are going to occur, and yet sufficient 
results may still be obtained because the upwind schemes are more forgiving in areas 
where central-difference schemes would require more points. A possible disadvantage 
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of upwind schemes is that since their dissipation is built in to their formulation, there 
generally is no mechanism to directly decrease the amount of dissipation introduced by 
the scheme for a given flow condition. Different upwind schemes can have different 
amounts of dissipation. 

Most upwind schemes have implicit time marching. So even though they may capture 
all of the desired features on a coarser grid than the central-difference computer codes, 
the central-difference computer codes can execute on a finer grid, which can provide 
sufficient resolution, in less time than the implicit, upwind computer code can on a 
coarser grid [71]. It was on this basis that an explicit upwind code was to be developed. 
The choice was made to use a modified Runge-Kutta explicit time integration method, 
very similar to that used in most central-difference computer codes. If modifying the 
coefficients used in the Runge-Kutta method provides sufficient damping for the general 
cases of interest, then the project will be considered successful. 

The first topic covered in the dissertation is the development of the governing equa- 
tions. Starting with the Navier-Stokes equations in dimensional Cartesian form. These 
equations are converted to nondimensional body-fitted coordinates, and are modified to 
a thin layer formulation, for all three coordinate directions. Turbulent equations are 
developed by executing the Favr6 density averaging on the Cartesian form of the full 
Navier-Stokes equations. The Favr<§ density averaging was chosen because the governing 
equations are for compressible flow. Algebraic turbulence modelling was chosen to re- 
solve the turbulent flows. The turbulence equations were written in thin-layer body-fitted 
coordinates. The Baldwin-Lomax [72] algebraic turbulence model was chosen to provide 
the eddy viscosity values used in the turbulence equations. An explanation of this model 
and the equations that are used is provided. 
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The two different types of upwind solvers are presented in the Numerical Schemes 
chapter. A brief explanation of van Leer’s [23] flux-vector splitting and Roe’s [25, 26] 
flux-difference splitting are presented. These two methods are dependent upon the method 
used to provide values for the splitting of the fluxes. The extrapolation method utilized in 
this work is the Monotone Upstream-centered Schemes for Conservative Laws (MUSCL) 
type differencing, which is subsequently explained. Then the time integration method 
used to advance the solution is provided, followed by a definition of the time step. At, 
and the different types of implicit residual smoothing, and implicit corrector smoothing. 

Chapters 2 and 3 provide the governing equations and the solution methods employed 
in this work. Chapter 4 explains multigrid acceleration and how it is employed in 
accelerating a numerical solution to a steady-state. It also provides a schematic of a 
computer listing to show an efficient method of programming this acceleration technique 
that allows grid independent subroutines. This is important because it sets the foundation 
for expanding the computer program to have multi-block capabilities. 

The multi-block structure developed for this computer program is also presented in 
Chapter 4. How it is implemented and a schematic of its incorporation into the computer 
listing is provided. The same chapter addresses the issue of the multi-block multigrid 
interactions, as well as provides the sequence for one multigrid cycle incorporating 
smoothers and multiple blocks. 

The boundary conditions used for this computer program are explained in Chapter 5. 
Along with the standard boundary conditions, the interface requirements are provided as 
well. The cases studied are provided in Chapter 6. They were chosen to validate the 
computer program and show its flexibility to handle different geometrical configurations. 
The final chapter contains concluding remarks about the cases studied and the success 
of the computer code. 
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Chapter 2 GOVERNING EQUATIONS 


2.1 Full Navier-Stokes Equations 


This work is directed toward solving compressible fluid flows for various configu- 
rations and physical cases. Many global properties of fluid flow can be obtained from 
simplified equation sets, but as the interests focus to areas closer to the body and global 
properties are not the only interest, more physics are needed in the governing equations. 
This is especially true where there are geometry changes in a body or in investigating 
wake regions and shear layers. It is these cases that require equations which include 
viscous phenomenon. For early researchers in CFD, the equations of choice were the 
potential equations coupled with boundary-layer equations. As computer equipment im- 
proved the more versatile Navier-Stokes equations became the governing equations of 
choice, and they are chosen as the governing equations of fluid motion for this work. 
Although they are computationally more demanding, they eliminate the approximations 
and restrictions implied in the older potential flow-boundary-layer approaches. 


The full Navier-Stokes Equations in dimensional Cartesian coordinate notation can 


be written as; 
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and 


The dimensional quantities, p, u, v, w, p, and E are the density, Cartesian velocities, 
pressure, and total energy per unit volume, respectively. The dimensional shear stress is 
represented by r, and q represents the dimensional heat flux. Invoking Stokes hypothesis, 
the shear stress terms become 
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where p is the dynamic viscosity. The formulation for the heat flux term is 
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where T is the temperature, and the coefficient of thermal conductivity, k, is given by 

l-ht ( 2 . 1 . 7 ) 

Pr 

where c p and P r are the specific heat at constant pressure and the Prandtl number, 
respectively. These equations are coupled with the perfect gas equation of state: 

p = pikf (2- 1 - 8 ) 
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where & is the gas constant. The assumption of a calorically perfect gas is made allowing 


c v T = e; 


( 2 . 1 . 9 ) 


where c v is the specific heat at constant volume, and e, is the internal energy. Thus, 


E = p 
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where the potential energy (PE) and other terms are assumed to be small or constant and 
have a negligible effect; therefore 
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and 7 = I s which is the ratio of specific heats. 

It is more convenient to have the equations in nondimensional form, because it 
allows for easier scaling of the flow case of interest, especially in investigating viscous 
flows. Doing so allows the independent variation of such parameters as the Mach 
number, Reynolds number, and Prandtl number, so that the computational results can 
be generalized and not be restricted to a specific geometrical configuration. Also, the 
flow variables are normalized so that their values will be between certain limits. Defining 
T, pref, T re f , o re /, and p Tef to be the reference length, density, temperature, speed ot 
sound, and dynamic viscosity, the following nondimensionalizations were employed; 
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These nondimens ionalizations allowed the governing equations to be rewritten as 
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and likewise for G, G„ H, and H v . Here H, M re/ , R ref , and 3 are the total 
Enthalpy per unit volume, reference Mach number, reference Reynolds number, and 

7-1, respectively. 

A Cartesian coordinate system may be ill-suited for many types of geometries. It 
can cause an inefficient use of points, and it can be very cumbersome to implement the 
proper boundaries for solving a given configuration. Using a curvilinear grid permits a 
better fitting grid around the body, and a more efficient use of cells. The cell sizes can 
change so that in regions of very small gradients, away from the body, large cells can be 
used. This allows more cells to be used near the body where the flow gradients are larger 
and require more cells for accuracy. The body-fitted coordinates can then be transformed 
into a computational domain that has equal sized cells. There are advantages to this, 
such as having the body surface selected as a boundary in the computational domain, 
allowing for easy application of the boundary conditions. Performing the transformation 
from Cartesian coordinates to curvilinear or body-fitted coordinates, f, rj, and ( is 
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accomplished by the following formulations: 
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It can be shown that the metrics are: 
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where the sub-characters x, y. and r on (, r/, and C represent partial derivatives of 
the body-fitted coordinates relative to the sub-characters. Here J is the Jacobian of 
transformation; 
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These formulations allow the governing equations to be written as; 
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where IK V , and U r are the contravariant velocities defined as, 
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Due to their complexity, G v . and H v , are provided in Appendix A. 

2.2 Thin-Layer Navier-Stokes Equations in Body-Fitted Coordinates 

Many fluid flow cases do not require the full Navier-Stokes equations. Generally there 
is a surface that has its normal perpendicular to the main stream wise flow. It is on this 
surface that a boundary layer will develop; therefore the cross flow and cross derivative 
viscous terms may be so small as to be considered negligible in comparison to the other 
flow terms. To sufficiently capture the boundary layer, many points are required near the 
boundary, thus allowing the governing equations to accurately predict the large gradients 
in the boundary layer. In some cases the physics of the fluid flow may require the full 
Navier-Stokes equations, but if there is not a dense packing of points in all directions the 
gradients of the flow will still not be captured, thus producing the same solution as that 
of the thin-layer equations. A compromise is to have the thin-layer approximation to the 
Navier-Stokes equations in all three coordinate directions. One reason is that it allows 
for generality in different flow cases, meaning that whichever direction is going to have 
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a boundary layer develop, the computer program is already capable of accommodating 
it. Plus, if there is significant cross flow and cross derivative flow characteristics, then 
the computer code is ready to accommodate them as well, provided there is sufficient 
packing of points in the required directions. 

By examining each viscous flux individually and eliminating the cross derivative 
terms in that flux, the three viscous, thin-layer Navier-Stokes fluxes in body-fitted 
nondimensional form are: 
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and e = j — jyp~' These three fluxes can replace the F v , G v , and H v in equation (2.1.18), 
which is the full Navier-Stokes equation in body-fitted nondimensional form. All of the 
other terms in that equation will remain the same. 


2.3 TVirbulence 

There are very few flows, outside of academic cases, that are laminar and attached. 
Most laminar cases are separated. The majority of flow cases are turbulent, which may 
have attached flow even though the laminar case would be separated. In other cases the 
turbulent flow may be separated as well. The larger the separation region, the greater the 
difficulty in resolving the flow. There are two ways turbulent flows can be resolved. One 
is by direct simulation, which uses the Navier-Stokes equations and directly solves for the 
different turbulent scales, down to the mean free path. This requires a tremendous number 
of points and has thus far been limited to a very small set of problems, for which the 
Reynolds numbers are in the range of one to three thousand [73]. The second approach is 
to use a turbulence model that will account for all of the different turbulent scales. This 
approach requires the use of the Reynolds averaged Navier-Stokes equations. Turbulence 
models can range from algebraic eddy viscosity models, which are the simplest, to second 
order closure models involving a minimum of seven additional differential equations 
which must be solved simultaneously with the original governing flow equations. 

Turbulence modeling was the approach chosen for the current work. The first step was 
to obtain the Reynolds averaged Navier-Stokes equations. Since the problems of interest 
are for compressible flows, the Favrd density averaging approach was selected [73, 74]. 
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The following formulations were used to perform the Favrd density averaging: 
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Using this formulation allows the continuity equation to be written as 
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Utilizing the average turbulent kinetic energy, A\ and the fluctuating kinetic energy, 
A, where 
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it can be shown that 
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This in turn yields 
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and is defined as the Reynolds stress tensor. Implementing the previous relation into the 
energy equation produces 
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The following terms need to be accommodated in order to solve the governing equations 
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The first term, pu,h, is generally considered the turbulent heat flux, and is often modelled 


as 
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where Jt T is the coefficient of turbulent thermal conductivity [73]. The second term, r* 
will be modelled using the Boussinesq s assumption [73, 75]. 


_R _ 


r“ = -puiUj = p T 


duj du, 2 du 

\_dx x dx } 3 dim u 


opA 8 X} 


(2.3.13) 


where pi is the turbulent eddy viscosity and is related to &t by 


fc T = 


CpPT 

PrT 


(2.3.14) 


R 

with P r T being the turbulent Prandtl number. It is common practice to combine r i} 
with Ti } to give 


rj, = + PT) 


duj du, 2 du m ^ 

[dx, dxj 3 dx m v 


- 3 TfKSij 


(2.3.15) 
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The remaining two terms and pu x K are either considered negligible and therefore 
accommodated by the turbulence model, or they are solved for as dependent variables in 
higher order schemes. For the present work these terms are assumed to be accommodated 
by the model. In executing the Favr£ density averaging, the equation of state needs to 
be included, which in turn becomes 


P - (7 - 1 ) 



p -^-~pK 


( 2 . 3 . 16 ) 


Thus, the turbulent full Navier-Stokes equations in dimensional Cartesian form become; 


dq d(f - fv) djg-jy) d{h - K) = Q 
dt dx dy dz 
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du 1 du 1 
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dv 1 du> 

51 + w 


hv = (fi + ht){ 


2\c.dw du _ dvl 1 

3 [^57 57 5y J + 

s(# + i) +:; (^ + t) + 

•»! [ 2 §? - if - fe] + “sJ+^T 


(/c+/ct) dT 


(2.3.20) 


the has been removed for clarity, and the lower case letters q , /, / t >, 9-, 9v, h, and h v 
are used to signify a difference between the other equation sets. 


2.4 Turbulence Modelling 

For the present work an algebraic turbulence model was chosen. Although these 
models generally do not resolve separated regions very well, they are easy to implement 
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and provide reasonable results for unseparated flows. Other turbulence models, such as 
the two equation models K - e, K - w, and K- r, are generally more applicable, 
especially in wake regions, but they require a certain amount of adjusting for a particular 
class of problems and are more complicated. The second order closure models are 
still in their early stages and are not commonly used for complicated geometries and 
configurations. It was decided that an algebraic turbulence model would be sufficient 
for the flows that were to be tested, and that it would provide a good starting point for 
validating the turbulence equations. More complicated models could be incorporated in 

the future. 

In using an algebraic turbulence model, the turbulent kinetic energy term. A, is 
dropped, because there is no mechanism in an algebraic turbulence model that can account 
for K. For consistency it is also dropped from the equation of state. Thus, writing the 
turbulent full Navier-Stokes equations in nondimensional Cartesian form produces; 


QQ d(F-F„) a( 6 -c„) 

dt + dx dy 



= 0 . 


( 2 . 4 . 1 ) 
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h = ^L(p + w){ 
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where 


a = 


2 ( JL . f2i\ 


(2.4.6) 


(3{fl + mt) ' ^rT / 

In comparing with the thin-layer Navier-Stokes equations in nondimensional Cartesian 

form, the only difference is between the viscous fluxes, F v , G v , and H v , where the 
viscoscity, p, is relaced by the sum p+pi, and jjr is replaced with /J(;I ^ T) (£ + #)• 
Therefore, to obtain the turbulent thin-layer Navier-Stokes equations in body-fitted 
nondimensional form, these two simple changes need to be made to the laminar, thin-layer 
Navier-Stokes body-fitted nondimensional equations. 


2.5 Baldwin-Lomax Algebraic Turbulence Model 

The Baldwin-Lomax algebraic turbulence model [72] is a two layer eddy viscosity 
model. The inner layer eddy viscosity model is the Prandtl-van Driest formulation defined 


as: 


MT inner 


pi 2 M 


H-ref 

M re f 


where 


l = kiy 



| w | = the magnitude of the vorticity. 



(2.5.1) 


(2.5.2) 


(2.5.3) 


and 

(2.5.4) 

Mref 

where, y is the distance normal to the wall, k\ is von Karman s constant (0.4), A 
is taken as a constant (26.0), and io max is the maximum vorticity along the coordinate 
direction normal to the wall. For a particular wall location, p w and p w are the values of 
density and molecular viscosity at the wall, respectively. The original Baldwin-Lomax 


y = y 


Pw^max 
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algebraic turbulence model did not use uw, but instead suggested using the shear stress 
at the wall, r w . This can produce erroneous values for the turbulence model if the 
flow is separated; therefore using u ma x has been found to be more reliable, and if there 
isn’t any flow separation on the wall, uw can be shown to approximately equal r„. 
The second part or outer layer of the Baldwin-Lomax algebraic turbulence model is the 
Clauser Formulation, given as; 


outer — 


A 2 Cep P F wake-^K LEB 


Rref 

M r ef 


(2.5.5) 


where 


( 


Vruax^n 


wake 


or 


Cwk y-max^dt f 

F max 


the smallest of 
the two values 


(2.5.6) 


with the Clauser constant K 2 = 0.0168, C cp = 1.6, C wh = TO for transonic flow, and 


V dif = (v / u 2 + w 2 + ^ 2 ) | 



+ V 2 + w 2 



(2.5.7) 


along the coordinate perpendicular to the surface at a particular wall location. For example 
the difference is along a constant x-location, if x is the streamwise direction. The value 
j /max corresponds to when F(y) = F{y) max , where 


F{y) = y|w| 



(2.5.8) 


In the wake region the exponential term for the previous equation is set equal to zero 


The Klebanoff intermittency factor is given as; 


Fkleb 


1+5.5 

y Ckleb 

6 “ 

Dmax 




(2.5.9) 


where Ckleb = 0.3. 
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To obtain the eddy viscosity, both the outer and the inner formulations are computed 
for that particular streamwise coordinate location. Then a comparison is made between 
the two viscosities to determine the location, starting from the wall or slip line, where the 
inner eddy viscosity value becomes larger than the outer eddy viscosity value. It is at this 
location that one switches from using the inner eddy viscosity model s values to using the 
values from the outer eddy viscosity model. The final eddy viscosity values are used in 
conjunction with the dynamic viscosity when the viscous flux derivatives are computed. 

2.6 Euler Equations 

The Euler equations are obtained by the simple elimination of the viscous terms, 
F v , G v< and H v , from the nondimensional, body-fitted, thin-layer Navier-Stokes equa- 
tions. This is easily accommodated in the computer code by having the viscous fluxes 
evaluated in a separate subroutine. 
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Chapter 3 NUMERICAL SCHEMES 

Upwind solvers were chosen to determine the inviscid fluxes. They gather informa- 
tion from both sides of an interface and then, based on the characteristic directions, a 
blending of the gathered information is performed. In this process if there is supersonic 
flow passing through a cell face then information from only the upstream side of the cell 
face is used. If the flow is subsonic, information from both upstream and downstream 
of the cell face is used. This process is based on a one-dimensional analysis, and thus 
assumes that information passes normal to the cell face. There have been studies to 
use true multi-dimensional characteristic directions, but they have met with only limited 
success [76]. The general approach is to perform one-dimensional analyses in the three 
coordinate directions independently, and then sum the results from the three directions. 
This was the method adopted for the present work. 

The viscous fluxes at the cell interfaces are evaluated using central-difference oper- 
ators, because these operators are better suited for these terms, especially evaluating the 
second order derivatives. 

3.1 van Leer’s Flux-Vector Splitting 

Flux-vector splitting is similar to the method of characteristics because it attempts to 
establish zones of influence and dependence in the flow field. Each flux-vector is split 
into a forward flux-vector and a backward flux-vector, allowing upwind differencing to 
be used for the spatial derivatives of the split fluxes. In 1982, van Leer [23] introduced 
a flux-vector splitting scheme that was designed to meet seven requirements. These 
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requirements were to provide, among other qualities, continuous split flux-vectors, so 
there would be smooth transitions when eigenvalues changed signs, and stationary shock 
structures within no more than two cells. These are qualities that the previous flux-vector 
splitting techniques lacked [20]. 

Following Ref. [23], the flux vectors F, G, and H can each be split into two vectors, 
a forward flux-vector based on non-negative eigenvalues, and a backward flux-vector 
based on non-positive eigenvalues. 

F = F + + F-, G = G + + <T, H = H+ + H~ (3.1.1) 


For local supersonic Mach numbers: 

Mi> 1.0, F? = F h Ff = 0 

Mi < —1.0, F, + = 0, Ff = Fi 
where / = f r/, and ( to indicate the three coordinate directions. 


( 3 . 1 . 2 ) 


For subsonic local Mach numbers, IM,I < 1.0 (in general notation for body-fitted 
coordinates [77]), a local scaled contravariant velocity component, u t , is defined as 

t T u + l,v + l'W l = (rl( (3.1.3) 

+ '? + '! 

where the local Mach number is given as 


Ui 

Mi = — 
a 


and a is the local speed of sound. The fluxes are: 
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where. 


and 


Here, 


ln = 


In 


ft + ', + '» 


n = x, y, z 


(3.1.6) 


fL„ = ± !) 2 


(3.1.7) 


' -0u] ± 2/3u/g + 2a 2 ( u 2 + + u.’ 2 

- + 2 
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± 

energy 


( 3 . 1 . 8 ) 


Fc = F F t) = G F< = H 


(3.1.9) 


= It U, = V U(_= w 


( 3 . 1 . 10 ) 


and 

0 = 7 -l (3.1.11) 

The “+” indicates the forward flux and the indicates the backward flux. 

Carrying out a central difference on each flux vector at the cell center gives: 

RHS = 


-Ad F+-F+ + F7 l -F 


+ G+ +i -G;_, + G- +i -G-_ i 

+ HU - H U + V ” H *-i 1 


(3.1.12) 


The present formulation, when applied to transonic and low supersonic flows, does 
not require the use of flux limiters for essentially oscillation free shocks. This was noticed 
by Anderson, Thomas, and van Leer [78], von Lavante and Haertl [79], Melson and von 
Lavante [80], and Cannizzaro, von Lavante, and Melson [81] and was explained in more 
detail by van Leer [23]. 
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3.2 Roe’s Flux-Difference Splitting 


Roe’s flux-difference splitting is an upwind scheme that approximates the Riemann 
problem at an interface between two cells by Roe’s averaging procedure [82]. The 
spatial derivatives across a cell interface (for example, the ^-direction) can be written 
conservatively as a flux balance across the cell in the form. 

* 2 . -- F j+ , - F i ( 3 - 2 -D 

*+T * 2 

where F is the numerical flux-vector at the corresponding cell interface. FoUowing Ref. 
[83], the cell interface flux is evaluated as 


F, +i = \\Fr{Qr} + F L {Q L ) - |/i|(A<?)] (3-2-2) 

Fl and Fr are the flux vectors computed from the left and right states, and A is the 
Roe averaged flux Jacobian matrix 



where 



(3.2.4) 


and 

A\ = Ss\A\Sl 1 ( 3 . 2 . 5 ) 

The A refers to the difference between the state variables on the left and right sides of 
the cell interface (For example, AQ = Qr - Ql)- and S" 1 are the left and right 
eigenvector matrices, respectively, for the direction, and A is the diagonal eigenvalue 
matrix. The present notation was adopted from Ref. [83], because it provides a simple 
programming strategy of these expressions. The ~ ’s refer to Roe averages computed as 

. _ UL\fPL_ + (3 2 . 6 ) 


34 



The last term in equation (3.2.2), 


(Qr-Ql) is a damping term due to the upwind 


character of the scheme and is given in detail in Ref. [82] as 


where, 
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(3.2.7) 


( 3 . 2 . 8 ) 


(3.2.9) 


with ^ = yM + / 2 + ^ , for / = ^ p, or C. and h representing the Roe averaged 
enthalpy. 


3.3 MUSCL Type Differencing 

Rather than determine the values of the inviscid fluxes at the cell centers and 
then extrapolate them to the cell interfaces. Monotone Upstream-centered Schemes for 
Conservative Laws (MUSCL) type differencing was used to determine the flux values at 
the cell interfaces. The dependent variables are extrapolated to the cell interfaces and 
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from these extrapolated values the fluxes at the cell interfaces are computed. This can 
be seen in Fig. 3.3.1. Thus, 

i) = 

a U= G+ ( Q k) G ^ = G "(<W> 03A) 

Where the “+” and are used to distinguish between the two directions of extrapolation. 
This approach provides a more accurate blending of the fluxes, because the blending will 
be based on the flow values at the interface rather than on some type of weighted averaging 
of the flow values from the separate cell centers. 

In many cases the primitive, or non-conservative, variables are extrapolated in the 
MUSCL approach rather than the conservative variables. 


i + 1/2 



Figure 3.3.1 Schematic of MUSCL Type Differencing 
The extrapolation procedure was written in the k - scheme formulation, where 

q t+ 1 = q, + ^^[(1 - K)V t + (1 + k) A<] (3.3.2) 

with V, = q t - qi-i and A t = 9l+1 - q x . The value of k determines the spatial accuracy 
of the extrapolation; k = -1 is pure second order accurate upwind, k = 0 is Fromm s 
(1968) scheme [84], which is second order accurate upwind biased, k = 1/3 is third 
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order accurate, upwind biased (it is less than third order accurate for multidimensional 
computations), and k = 1 is the second order accurate central difference. If “ swtch ” is 
set equal to zero, then the extrapolation is first order. 

3.4 Time Integration Method 

There are two types of time integration schemes used to advance numerical calcula- 
tions to steady state solutions. One is implicit and the other is explicit time integration. 
Implicit time integration schemes [11-18] require more computational work per time 
step, but they allow the time steps to be large and aid in the propagation of boundary 
information because of their elliptic nature. Many implicit schemes involve the solu- 
tion of either scalar or block tri-diagonal matrices, or they approximate these matrices 
with bi-diagonal or even diagonal matrices [85]. Explicit methods usually require at 
least two computation stages, such as MacCormack’s [8] predictor-corrector method, or 
they can have many stages, such as a Runge-Kutta multistage scheme. Modem imple- 
mentation of Runge-Kutta methods can be found in the work of Jameson, Schmidt, and 
Turkel, [32], Jameson and Baker [33], and Jameson [34]. More stages generally permit 
a higher CFL number, and better error smoothing properties. With each stage, the flux 
evaluations are computed. After completing all of the stages, the numerical solution is 
advanced one time step. Most upwind methods use an implicit time integration technique 
because it provides better smoothing and the fluxes are evaluated only twice per time 
step (once for time level “n”, and once for time level “n+1”). Thus, if a multistage 
method is employed, more stages may be required to provide comparable smoothing, 
causing more flux evaluations, which could lead to a higher computational effort than the 
implicit time integration approach. The explicit time integration schemes are generally 
central-difference schemes because the flux evaluations cost much less than the upwind 
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approach; therefore more stages can be executed at a cheaper computational cost than 
using an implicit time integration. 

The current work was directed at having an upwind explicit computer program. 
Upwind methods were chosen for their accuracy, and explicit time integration was 
chosen in an attempt to capitalize on what the central-difference explicit schemes have 
accomplished, in terms of computational speed. Also, in considering a multi-block 
computer program, there will be many flow configurations that will not allow the implicit 
sweeps to be conducted across all of the blocks simultaneously. Instead, each implicit 
sweep would have to be performed one block at a time, reducing the effectiveness of 
propagating boundary information across the entire domain. Although this will change the 
convergence rate, it has not been reported to cause a divergence of the numerical solution. 

The basic characteristics that are desired for the explicit time integration are good 
high frequency error damping, a minimum amount of computational effort, and robust- 
ness. The damping qualities are required for two reasons. One obviously is because the 
computer program will converge properly and expediently. Second, if good high fre- 
quency error damping is not achieved, multigrid acceleration will not perform properly. 
The amount of computational effort required is important, because if too many stages are 
required, it would cost less to use implicit time integration. If the scheme is not robust, 
it will not be generally applicable to various flow configurations, which would defeat the 
main purpose of developing this computer program. 

To obtain these characteristics it was determined that the ability to adjust the stability 
range and amplification factor would provide avenues by which the explicit technique 
could be tuned to satisfy the necessary damping qualities. The multistage Runge-Kutta 
method has been modified for many of the explicit central difference codes [32, 86], and 
therefore was chosen as the time integration technique for the proposed computer code. 
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The classical fourth order Runge-Kutta time integration scheme was used by Jameson, 
Schmidt, and Turkel [32] for solving a finite-volume formulation of the Euler equations 
(It provides fourth order temporal accuracy for linear equations and second order accuracy 
for non-linear equations). This approach requires storing the solution at each stage as 

follows; 

Q' =Q"~y RHS{Q n } 

=«--T RHS{Q ' ] ( 3 . 4 . 1 ) 

Q 3 = Q* - At RHS { Q 2 } 

q»+ 1 = Q" - J. [RHS{Q n } + 2 RHS{Q'} +2 RHS{Q 2 } + RHS{Q 3 }} 

This would be very memory intensive for a three-dimensional computer code (where 

RHS represents the flux evaluations). In 1983, Jameson and Baker [33] presented the 

following four stage Runge-Kutta scheme; 

Q 1 = Q n - aj A t RHS{Q n } 

Q 2 =Q n -a 2 AtRHS{Q 1 } 

* (3.4.2) 

Q 3 =Q n - ot 3 At RHS{Q 2 } 

Q n+i = Q”-a 4 A t RHS{Q 3 } 

The coefficients ai, and a 3 could be independently varied to obtain a range oi 
stability and amplification properties. The coefficient for the last stage, a 4 , must equal 
one for consistency, and also to provide at least first order accuracy in time. This new 
approach was also adopted by others [87], Initially, researchers still used the standard 

four stage scheme coefficients of aj = i, o 2 = ^ Q 3 = 5. and a 4 = which can be 
obtained from a Taylor series expansion. Attempts to obtain more desirable stability and 
amplification characteristics were pursued by varying the coefficients and the number of 
stages [42, 88, 89]. These approaches were soon investigated for upwind schemes by the 
present author [35] and others [36, 37]. The coefficients developed by coworkers [35] for 
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this computer program are presented for various extrapolations in the following tables 3.1 
3.2. The linear wave equation was used as the model equation. A full explanation of the 
approach and detailed results can be found in the Ph.D. dissertation of Alaa Elmiligui 
[90], who was a coworker in the development of the computer program. 


Table 3.1 Multistage Coefficients for First and Pure Second Order Schemes. 


IBMUHiW 

First Order 

Pure Second Order | 

Number of Stages 

Nui 

nber of St< 

rges 

Multistage Coefficients 

2 

3 

4 

2 

3 

4 

a l 

0.220 

0.105 

0.056 

0.220 

0.150 

0.091 

C*2 

1.000 

0.325 

0.152 

1.000 

0.400 

0.240 

03 


1.000 

0.340 


1.000 

0.420 

a 4 



1.000 



1.000 


Table 3.2 Multistage Coefficients for Fromm and k = 1/3 Schemes. 


1 •' : :*.• ‘i: 1 •: ; ”* * = • 

Fromm 

AC = 1/3 1 


Number of Stages 

Nur 

nber of Su 

lges 

Multistage Coefficients 

2 

3 

4 

2 

3 

4 

Oil 

0.420 

0.210 

0.110 

0.460 

0.220 

0.135 

a 2 

1.000 

0.440 

0.255 

1.000 

0.480 

0.260 

Ot 3 


1.000 

0.46 


1.000 

0.440 

a 4 



1.000 



1.000 


3.5 Local Time Stepping 


Local time stepping allows each cell to advance in time at its own or local stability 
limit. This approach provides for faster signal propagation, which produces faster 
convergence to a steady state solution. The local time step At is based on the Courant- 
Friedrichs-Lewy (CFL) stability limit. It is calculated as follows; 


1 1 1 1 1 
~At - ~KTt + At v + At c + At 1 ’ 


+ 


At* 


+ 


1 

At* 


(3.5.1) 


40 




















































where 


—rj— > A/ = |(//| + (3.5.2) 

At/ y 

with Ui = u/j. + u/y + w/ z being the contravariant velocity for the “/” direction, a the local 
speed of sound, A / the eigenvalue, and 4>i = yjp x + fy + The viscous contributions 
to the time step are: 


Aff “ ' P 


<t>i 


|max^-, | + jjdWsI + + IVd) 


(3.5.3) 


The first three terms on the right hand side of equation (3.5.1) are due to the stability 
limitation on the inviscid flux, while the last three terms are due to viscous flux stability 
limitations. The viscous time step limitation terms, make the scheme more robust 
on fine viscous grids for boundary-layer type flows [91, 92]. 


3.6 Implicit Variable Coefficient Residual Smoothing 

The purpose of residual smoothing is to reduce the magnitude of any spikes in the 
residuals. The residuals are generated by executing the flux differentiations, and are 
used to adjust the values of the dependent variables. These adjustments are necessary 
for the dependent variables to obtain their correct values for the given flow conditions. 
Residual smoothing can reduce high frequency errors, which is beneficial for multigrid 
acceleration techniques. Reducing these errors allows the restriction process in multigrid 
to be implemented without causing aliasing (Aliasing is further explained in the chapter 
containing multigrid). Implicit residual smoothing can increase the stability region and 
enhance the damping properties of a multistage time-stepping scheme. The formulation 
for three-dimensional problems is usually applied in the form, 

(/ - 0*V*A*)(/ - /3„V,A„)(7 - (3 C V C A C )R* = R (3-6.1) 

where R is the residual, and V and A are the standard backward and forward difference 
operators relative to the ?/. and ( directions. The coefficients 0 ^ , dii- and 0^ can be 
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constant, which applies the same magnitude of damping on every residual value. That 
approach is useful but not as general as having variable coefficients. Variable coeffi- 
cients self-adjust to produce damping only where needed; therefore variable-coefficient, 
residual-smoothing can provide better convergence than the constant coefficients without 
tuning. The two-dimensional variable coefficient method presented by Swanson, Turkel 
and White [93] was extended to three-dimensions as follows: 


A 



£IL * YW1 

CFL* A* + A„ + A c )\ J 


(3.6.2) 


where / = £, rj , and (• This formulation makes A a function of the grid aspect ratio and 
the spectral radii. A,. The ratio of is the CFL of the smoothed scheme to that of 
the basic explicit scheme CFL*. The optimum ratio was found to be §f^ » 2.0 [93]. 
This operator was applied before each Runge-Kutta stage. 


3.7 Implicit Corrector Smoothing 

Corrector smoothing was applied using constant coefficients. This technique was to 
provide a better multigrid correction value from the coarser meshes, by eliminating any 
erroneous spikes in the correction data. 
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Chapter 4 MULTIGRID MULTI-BLOCK 
AND THEIR INTERACTION 

4.1 Multigrid 

Multigrid is a technique used to accelerate the rate of convergence. It has been 
adapted to solve the ordinary and partial differential equations found in fluid mechanics. 
Acceleration is achieved by attacking the low frequency errors, which generally are not 
well damped by the equation solver. Most equation solver techniques can damp the high 
frequency errors, but require more iterations to damp out the low frequency errors, and 
as the number of mesh or grid points increases, it takes the original equation solver many 
more iterations to reduce the low frequency errors. The number of iterations increases 
non-linearly with the number of cells. True multigrid performance does not diminish 
with an increase in the number of grid points. Hence, multigrid can provide a converged 
solution in the same number of cycles as a grid that contained only every other point. 

The multigrid technique reduces the low frequency errors by solving a set of 
governing equations on successively coarser grids. Thus, what was a low frequency 
on a fine grid becomes a higher frequency on a coarser grid. In multigrid a fine grid 
that has had every other point eliminated in all directions is defined as a coarser grid. 
Thus, the equation solver is again working on high frequency error, for which it is best 
suited. Information gained from the coarser grids is used to reduce the low frequency 
errors on the finer grids. Generally three or four grid levels are used for a calculation. 
Also, running calculations on the coarser grids requires less computational work because 
of reduction in the number of points. 
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Multigrid acceleration techniques were originally applied to linear elliptic equations. 
They have since been modified to handle non-linear hyperbolic equations. Excellent 
developments of multigrid techniques can be found in references [94-96]. Multignd 
performs a certain amount of averaging of flow variables in restricting the values from a 
fine grid to a coarser grid, and in prolongating the correction to the flow variables from a 
coarse grid to a finer grid. Thus, it is better suited for elliptic flows than for parabolic or 
hyperbolic flows, due to the way information is physically propagated in the flow field. 
At the present, multigrid techniques have had their biggest impact in transonic flows, 
but modifications are being introduced to handle other demanding flow cases, such as 
hypersonic flows. 

A brief explanation of multigrid will be presented in two sections. The first section 
will explain multigrid methods for linear equations. This will provide the basis for the 
second section, which will explain multigrid methods with non-linear equations. 

4.1.1 Linear Equations 

Consider the problem 

l h lj h — f h (4.1.1) 

where L h is a linear, finite-difference operator on a grid, g h , and h is the cell spacing. 

The forcing function, /\ is known and U h is the solution to the problem on the grid 

with spacing h. If we take u h as an approximation to U h with an error of V h , i.e. 

V h = U h - u h (4.1.2) 

then equation (4.1.1) can be written as 

L h (u h + V h ) = f h < 4 - L3) 
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Since L h is a linear operator, this can be written as 

L h (u h ) + L h (v h ) = f h ( 4 - 1 - 4 ) 

If V h is a smooth function, meaning it does not have any high frequency errors, it can be 
represented on a coarser grid, g 2h , with spacing 2 h, which has twice the spacing between 
points as the grid with spacing h. The grid g 2h , is formed by removing every other point 
in grid g h . Therefore, g 2h € g h ■ Points are eliminated from g 2h to form g 4h and so 
forth to form g 8h , g l6h , etc. Each subsequent grid is a subset of the previous grid. (If a 
function is not smooth, aliasing will occur during the transformation of information from 
the finer to the coarser grids, thus preventing an accurate representation of data from the 
finer grid on the coarser grid.) 

It is possible to solve for an approximation to V h on grid g 2h , using the equation 
L 2hfj2hyh ^ = j2 h^fh _ L h u h ^ (4.1.5) 

where lf h is the restriction operator which transfers the values of a function from the 
fine grid to the coarse grid (An explanation of the restriction operator and how it is 
implemented can be found in Appendix B). If the coarse grid forcing function is defined 

as 

f 2h = ll h (f h - L h u h ^j ( 4 - 1 - 6 ) 

and the coarse grid error is taken to be 

y2h _ j2hy k ( 4 . 1 . 7 ) 

then 

L 2hy2h _ j2h ( 4 . 1 . 8 ) 

Since equation ( 4 . 1 . 8 ) is for a grid that is coarser than equation ( 4 . 1 . 1 ), the numerical 
evaluation of V 2h is much cheaper than the evaluation of V h on the fine grid. Once V 2h 
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is obtained, it is used to correct the fine grid iterative solution, u h , using 

(u h ) = (u h ) + I% h V 2h (4.1.9) 

The coarse grid to fine grid transfer operator, I% h , is the prolongation operator (An 
explanation of the prolongation operator and how it is implemented can be found in 
Appendix B). 

Since the form of equation (4.1.8) is the same as equation (4.1.1), it is obvious that a 
grid with spacing 4 h can be used to find corrections to the “solution” of the problem on 
the grid with spacing 2 h. Successively coarser grids may be used until a grid is reached 
which is so coarse that a direct solution may be used (or a nearly exact solution with 
only a few relaxation sweeps). The correction from the coarsest grid is then used to 
correct the correction on the next finer grid; and this is continued through successively 
finer grids until the finest level is reached and the approximate solution is updated. 

The usefulness of corrections obtained on a coarser grid is dependent on the smooth- 
ness of the fine grid error passed to the coarse grid. Hence, it is absolutely necessary that 
the high-frequency components of the error on the fine grid are reduced, if not completely 
eliminated. It is the responsibility of the smoother (usually a relaxation algorithm) to 
damp the high frequency components of the error. The removal of the low-frequency 
components of the error is unimportant for all but the coarsest grid since these frequencies 
can be resolved on the coarser grids where they become high-frequencies. If the high- 
frequencies are not damped, then the restriction operator will pass aliased information to 
the coarser grid and the entire multigrid scheme will cease to converge [97]. Obviously, 
the choice of smoother is critical to multigrid functioning properly. 
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4.1.2 Non-linear Equations 


The previous development of the multigrid scheme was for linear operators. Unfor- 
tunately, many problems in engineering are described by non-linear equations or sets of 
equations. This is particularly true in the field of Computational Fluid Dynamics (CFD). 
Because of the non-linear nature of the equations, the Full Approximation Storage (FAS) 
multigrid scheme [40] must be used (FAS is applicable to both linear and non-linear 
problems). A brief description of FAS follows and relies heavily on the description of 
multigrid for linear problems given in the previous section. 

In the development of multigrid for linear problems, the linearity of the operator 
was used to split the error out from the approximate solution as shown in the step from 
equation (4.1.3) to equation (4.1.4) above. If the operator is non-linear, this splitting is not 
valid. Instead, the derivation proceeds as follows, starting with the non-linear problem: 


jji jjh __ jh ( 4 * 1 . 10 ) 

Again, the substitution for the exact solution is made to give: 

£*(„* + V h ) = /‘ (4.1.1 1) 

Now, L h u h is subtracted from both sides of equation (4.1.11) to give: 

L h (u h + V h ) - L h (u h ) = f h - L h (u h ) (4.1.12) 

On the coarse grid, equation (4.1.12) becomes: 

L 2h (l 2 h h u h + V 2h ) - L 2h (l 2 h h u h ) = /?(/* - L h u h ) (4.1.13) 


As in the linear case, it is assumed that the error, V h , can be represented accurately on 
the next coarser grid as V 2h . If the second term on the left-hand side is moved to the 
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right-hand side, equation (4.1.13) can be written as: 

L 2h (u 21 ') = f a (4.1.14) 

where 

f U = jUtfl _ £ V) + £»(/»„») (4.1.15) 

and 

= /J* ti' 1 + V 2/l (4.1.16) 

The values of u 2h are obtained on the coarse grid and used to update the fine gnd 
solution using the following equation: 

Note that the prolongation term on the right-hand side of equation (4.1.17) is the correction 
to be applied to the fine grid solution. Examination of this term shows that the solution 
on the coarse grid is actually a solution to the originally posed problem and not just a 
correction to the fine grid solution as in the linear case. This is an important difference 
because it allows the use of the fine grid boundary conditions on all the coarse grids as 
well. As with the linear problem, the non-linear FAS scheme uses the same operator on 
all the grids. This of course simplifies the programming of the multigrid scheme. 

The following section describes the data structure and programming approach chosen 
in an attempt to efficiently code the multigrid acceleration technique. The structure was 
then utilized to allow the inclusion of multi-block flexibility with a minimum amount of 
changes to the existing computer program. Setting up the data structure in this manner 
allows considerable flexibility in how the subroutines are coded, as will be explained in 
the following sections. 
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4.1.3 Fortran Data Structure 


In explaining how the data structure was chosen, an understanding of how the memory 
is set up in a computer is needed. Despite the number of dimensions in an array, all 
data are stored in a one-dimensional array. During compilation, pointers are generated 
that indicate the register locations where different data can be acquired. This process 
is invisible to the user, but knowing how the memory process works allows one to 
exploit it to conserve memory space and write a general computer program that is easy 
to read. Knowledge of how the computer memory is arranged provides flexibility in 
how subroutines can be written, and be generic in terms of the different types of mesh 
topologies it can handle. 

It is a general practice in Fortran to use common blocks to dimension arrays, and 
transfer those arrays from the main program to the subroutines and from subroutine 
to subroutine. When using common blocks, the array dimensions are set in the main 
program and in all the subroutines at compilation, thus fixing the size ot the arrays in 
the subroutines. This becomes a disadvantage for a multigrid computer code, because 
of the numerous grid levels (generally a minimum of three). Each successively coarser 
grid level, for a three-dimensional computer code, has one-eighth as many grid points 
as the next finer grid level; therefore requiring much less memory space than for the 
finer grid level. The arrays could be made one-dimensional, and include all grid levels. 
Unfortunately this approach requires special counters and pointers to allow calculations to 
be preformed on the different grid levels. Another approach would be to have a separate 
three-dimensional array for each grid level, but it would be very cumbersome to manage. 
Furthermore, separate subroutines would be required to handle the different grid levels, 
which would obviously not be efficient. The most straightforward approach is to account 
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for all of the points that an array will need, including all grid levels, and set this array up 
as a one-dimensional array in a common block in the main program. However, instead of 
using common blocks in the subroutines, the main program should pass the information 
required by each subroutine through its argument list. When an array is passed through 
a subroutine argument, what is actually passed is the starting location, or address, where 
that array data can be found. Therefore, by passing to the subroutine the starting location 
and the dimensions of the array, for that particular grid level, the array can then be 
dimensioned in the subroutine during execution time, allowing the size of the array to 
change depending on the dimensions passed to the subroutine. The subroutine then deals 
with only that section of the array that the prescribed dimensions allow it to access. 

In the main program, the arrays used for multigrid are stored as one-dimensional 
arrays. An integer is then created whose values are the starting addresses for the multigrid 
data stored in the array. Generally, the fine grid data are stored first, then the intermediate 
grids, and then finally the coarse grid data. A schematic drawing of this storage sequence 

is shown in Fig. 4. 1.3.1. 


Starting address tor grid 8h 
Starting address for grid 2h 


I 



\ 

Starting address for grid h 


Starting address for grid 4h 


Figure 4. 1.3.1 Multigrid Storage Arrangement for Arrays. 

The starting locations and dimensions of the data arrays for each of the multiple 
grids are pre-calculated and stored in integer arrays as a function of the grid. A Fortran 
example is given below for the grid arrays containing the x, y, and z coordinates of the 
points in a three-dimensional computational domain. Several notational conventions are 
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used throughout the paper. The variable ngrid always refers to the maximum number of 
multigrid levels and grid level 1 refers to the finest grid. 


program multigrid 
parameter (ngrid=4) 


parameter 

parameter 

parameter 

parameter 

parameter 


j 1=33 , kl=17 ) 

/2, j2= (jl + l) /2 , k2=(kl + l)/ 
12 , j 3= ( j2 + l) /2 , k3= (k2+l) / 
( i4= (i3+l ) /2 , j 4= ( j 3 + 1 ) /2 , k4=(k3 + l)/ 

( i jkmax=il* j l*kl + i2*j2*k2 + 

i3* j 3*k3 + i4* j 4*k4 ) 


( il=97 , 
(i2= (il+1 
(±3= (i2+l 


2 ) 

2 ) 

2 ) 


common /mg/ istart (ngrid) , idim(ngrid) 
common /coord/ x ( i j kmax) y ( i jkmax) , z ( i jkmax) 


imax(l) = il 
jmax(l) = jl 
kmax(l) = kl 
istart ( 1 ) =1 

do 100 igrid=2 , ngrid 

imax(igrid) = ( imax ( igrid-1 ) + 1) / 2 

jmax(igrid) = ( jmax(igrid-l) + 1) / 2 

kmax(igrid) = (kmax ( igrid-1 ) +1) / 2 

istart ( igrid) = istart ( igrid-1) + 

imax (igrid-1) *jmax ( igrid-1 ) * 
kmax (igrid-1) 

100 continue 


do 200 igrid=l, ngrid 

call metrics (x( istart (igrid) ) ,y (istart (igrid) ) , 
z ( istart ( igrid) ) , . . - 
imax ( igrid) , jmax( igrid) , 
kmax ( igrid) , . . . ) 

200 continue 
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subroutine metrics (x,y, z, . . . , imax, jmax, kmax, . . . ) 
dimension x ( imax , j max , kmax ) ,y (imax, jmax, kmax) , 
z ( imax , j max , kmax ) 


4.1.4 V- and W -Cycles 

One multigrid cycle is started by performing work, or iterations, on the finest grid 
level, restricting that information to the next coarser grid, performing iterations on that 
grid level and then continuing on to a coarser grid level. Once iterations have been 
performed on the coarsest grid level, the correction information is prolonged to the 
next finer grid level. This continues until the last iterations are performed on the finest 
grid level, thus completing one multigrid cycle. The cycles are repeated until sufficient 
convergence is obtained on the finest grid. In this section, fixed cycles known as V- and 
W-Cycles will be described. 


In the present work, Fortran IF statements were avoided as much as possible. This 
led to a method of coding the multigrid cycles that relied heavily on DO loops. Basically, 
a standard V-Cycle can be broken into halves. The first half is the restriction part of the 
cycle going from the fine grid through the coarser grids down to the coarsest grid. The 
second half is the prolongation part of the cycle going from the coarsest grid up to the 
finest grid. An example is shown in Fig. 4. 1.4.2 for a four level multigrid. The circles 
indicate when iterations are performed on the given grid level, and the lines between grid 
levels indicate either a restriction or prolongation operation. Notice that the circle for the 
fine grid at the beginning of the cycle is omitted since the iterations on the fine grid are 
performed at the end of the prolongation section. This ensures that the last operations 
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Figure 4. 1.4.2 Schematic of Computation Sequence for a Four-Grid V-Cycle. 

in a multigrid cycle include updates on the fine grid. The control of the grid level is 
handled by DO loops as shown in the following, where ncycle is the total number of 
multigrid cycles to be performed. The fine grid is 1, the coarsest grid is ngrid\ iterate 
is the iterative solver, restrict performs the restriction operation, and prolong performs 
the prolongation operation. 


igrid=l 

call iterate ( . . . , igrid, . . . ) 


do 5000 icycle=l, ncycle 


do 1000 igrid=2 , ngrid, 1 
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call restrict ( . . . , igrid-1 , igrid, . . . . ) 
call iterate (..., igrid, ... ) 


1000 continue 


do 2000 igrid=ngrid-l, 1, -1 


call prolong ( . . . , igrid+1 , igrid, . . . ) 
call iterate ( . . . , igrid, . . . ) 


2000 continue 
5000 continue 


It is often necessary to perform more than one iteration on a given grid level to get 
the required smoothness in the error for multigrid to perform correctly. For simplicity, 
this iteration loop has been left out of the section of code shown above. 


A W-Cycle can be thought of as consisting of several components which are similar 
to V-Cycles but with varying ‘coarsest’ and ‘finest’ grid levels. This idea is shown in 
Fig. 4. 1.4.3, where a W-Cycle is expanded graphically to show its ‘legs’. This requires 
a simple coding modification to the V-Cycle program to allow W-Cycles. Another DO 
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Fine 




loop is added to perform the various legs of the W-Cycle and then the beginning and 
ending grid levels of the legs are varied to create the W-Cycle. For the four-gnd-level 
W-Cycle shown in Fig. 4. 1.4.3, the fine and coarse grids, as a function of the W-Cycle 
‘leg’, are shown in the following table. This W-Cycle has only four legs; the fifth leg 


Table 4.1. - Legs for Four Level W-Cycle. 


LEG 

FINE GRID 

COARSE GRID 

Heg 

iline(ileg) 

icoarse(ileg) 

1 

1 

4 

2 

3 

4 

3 

2 

4 

4 

3 

4 

(5) 

(1) 



is actually the first leg of the next W-Cycle, but it is shown here for completeness. A 
representative listing of the Fortran coding for the W-Cycle is shown below. The variable 
nleg is the number of legs in the cycle; for example, nleg = 4 for the W-Cycle shown in 
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Fig. 4. 1.4. 3. The variable ileg specifies the current leg of the cycle and varies from 1 to 
nleg. The bold characters show the changes from the V-Cycle code to the W-Cycle code. 


igrid -ifine (1) 

call iterate ( . . . , igrid, . . . ) 


do 5000 icycle=l , ncycle 

do 3000 ileg=l, nleg 


do 1000 igrid-i fine (ileg) +1, icoarae (ileg) , 1 


call restrict ( . . . , igrid-1, igrid, . . . . ) 
call iterate ( . . . , igrid, . . . ) 


1000 continue 


do 2000 


igrid=i coarse (ileg) -It ifine ( ileg+1 ) 


-1 


call prolong ( . . . , igrid+1 , igrid, . . . ) 
call iterate ( . . . , igrid, . . . ) 
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2000 continue 

3000 continue 

5000 continue 


One method often used to increase the rate of convergence when utilizing the multi- 
grid process is to execute what is called full multigrid (FMG). If the grid configuration 
contained enough cells to support a four level V-Cycle, FMG would start by just cycling 
through the coarse grids; therefore the initial V-Cycle may be only a two level V-Cycle, 
where the two grids used are the two coarsest grids of the possible four grids. Of the 
two coarsest grids, the finer is treated as a solution grid; therefore its multigrid forcing 
function is zero. After a sufficient drop in the residual the next finer grid will be in- 
cluded into the process, making either a three level V-Cycle or the coarsest gnd could 
be dropped from the cycle so that only a two level V-Cycle will still be used. It is the 
computer code operator’s decision whether a two or three level V-Cycle is used. This 
process is continued until the finest grid is incorporated into the V-Cycle. The objective 
is that convergence will be enhanced because on just the coarser grids the solution will 
set up faster and also require less CPU time; therefore, once the fine grid is included 
into the V-Cycle, the large flow field characteristics should be developed. The FMG for 
a V-Cycle is shown in Fig. 4. 1.4.4 and for a W-Cycle in Fig. 4. 1.4.5. 
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Figure 4. 1.4.5 Schematic of Full Multigrid Four Grid Level W-Cycle 


4.2 Multi-block Structure 

The multi-block concept is best explained by considering complex grid configurations. 

If one considers the case of an engine mounted on a wing/body configuration, or an 
intemal/external afterbody, or any other designs containing multiple solid surfaces, it 
becomes apparent that these designs cannot be handled with a simple single-block 
configuration. But, if each element of the configuration can be accommodated separately, 
the grid generation task becomes much easier. Plus, in some cases more than one type of 
grid topology may be required. One task is the grid generation of these configurations, 
and deciding what grid topology to employ. One body part may best be fit with an H-H 
topology, while another may be best fit with an H-O topology. A multi-block code which 
can handle multiple solid boundaries for various topologies and configurations, such as 
those mentioned in the first chapter, without requiring changes to the source code, allows 
one to grid a given configuration in a manner that provides an optimum topology for the 
different components of the grid. 

A multi-block computer code allows a complex configuration to be divided into 
sections or blocks. Each block can be defined as a six sided volume, and can be of 
a different grid topology. Each side or face of the block can have different boundary 
conditions than the other faces. Plus, each face can be divided into multiple segments or 
patches, with each patch having a different boundary condition. The blocks communicate 
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with each other through interface conditions. Interface conditions can also be used to 
allow a block to communicate with itself, such as across a wake cut for a wing calculation. 
The simplest interface condition is for the grid lines to have C 1 continuity at the interfaces, 
meaning the grid lines have to match one for one (be homogeneous) across the interface, 
and that the grid be continuous across the interface. More complex interface conditions 
would allow more grid points on one block face of the interface than on the other block 
face of the interface. The type of interface conditions that can be allowed depends on the 
sophistication of the interface routine. Note that the only way to maintain higher order 
accuracy and flux conservation across an interface is by having a continuous grid across 
the block interfaces. Further discussion on interfaces will be provided in the boundary 
condition chapter. 

The present multi-block code was designed on the premise that if it can accommodate 
multiple blocks, and each block can have various boundary conditions, then the source 
code can handle any geometrical configuration that can be represented with six-sided 
blocks. The block and boundary information is provided to the executable through 
an input file. This allows the source code to accommodate the different geometrical 
configurations without being altered. 

Once the code was developed to handle the multigrid format, the amount of work 
needed to include multi-block flexibility was greatly reduced; since the subroutines were 
grid level independent for the multigrid structure, they will be block independent as well. 
The only constraints on some of the subroutines will be the boundary condition ranges, 
which will be passed through the argument list. The last decision before making changes 
which would provide the multi-block capabilities, is whether to execute all multigrid 
levels in each block and then proceed to the next block (this method is called executing 
multigrid inside of multi-block), or execute the flux evaluations for all the blocks at a 
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particular multigrid level then proceed to the next multigrid level for all blocks (this 
method is called executing multi-block inside of multigrid). It was decided that better 
communication would be delivered if multi-block inside of multigrid was used. This will 
be explained further in the multi-block multigrid interaction section. 

4.2.1 Multi-Block Storage and Programing Strategy 

Basically, each of the subroutines in the multigrid computer code was designed to be 
independent of the multigrid level. This same characteristic makes each of the subroutines 
independent of the block passed to it in a multi-block environment; therefore, the memory 
allocation scheme described in figure 4. 1.3.1 is expanded to include multiple blocks as 
shown in figure 4.2. 1.1. 


Block 1 


Block 2 


Block 3 


Block 4 


Block 5 


Figure 4.2. 1.1 Multigrid Multi-block Storage Arrangement for Arrays. 
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To accommodate the extra memory needed for the multi-block strategy, an extra 
index is added to the pointers and parameters used for finding a specific location in an 
array and defining the size of the arrays, respectively. The changes to the previously 
shown multigrid code structure are indicated in bold print. 


c 

c 

c 


c 

c 

c 


c 

c 

c 


program multigrid 
parameter (ngrid=4) 

parameter (nblock=2) 

BLOCK 1 

parameter (ill=97, 
parameter (i21= (ill+1) /2, 
parameter (i31= (i21+l) /2, 
parameter ( i41= ( i31+l ) /2 , 
parameter ( i jkmaxl= (ill* j 

i3 1* j 3 1 


j 11=33 , kll=17 ) 

j 21= ( j 11+1 ) /2 , k21= (kll+1) /2 ) 
j31=(j21+l) /2, k31=(k21+l)/2) 
j41= ( j 31+1 ) /2, k41= (k31+l) /2 ) 
l*kll + i21*j21*k21 + 
k31 + i41* j41*k4i; ) 


BLOCK 2 


parameter (H2-65, j 12-49, kl2 33 ) 

parameter (i22=(H2+l) /2, j22=(jl2+l) /2, k22= (kl2+l) /2) 
parameter (132= (122+1) /2, j32=(j22+l) /2, k 32 = (k22+l) /2 ) 
parameter (142= (132+1) /2, j42= (j32+l) /2, k42= (k32+l) /2 ) 
parameter (ijkmax2= (H2*jl2*kl2 + i22*j22*k22 + 

i32*j32*k32 + 142*j42*k42 ) ) 


TOTALS 

parameter (Ijkmax = ijkmaxl+ijkmax2) 


common /mg/ istart (ngrid, nbloc) , idim(ngr id, nfcloc) 
common /coord/ x ( ij kmax) , y ( ijkmax) , z ( ijkmax) 


imax ( 1 / 1 ) = i 1 1 
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j max (1,1) = j 1 1 
kmax(l ,1) = kll 
istartd ,1) = 1 


do 100 igrid=2 , ngrid 

imax(igrid/l) = (imax ( igrid-1/ 1) + 1) / 2 

jmaxdgrid ,1) = ( jmax ( igrid-1, 1 ) + 1) / 2 

kmaxdgrid,l) = (kmax ( igrid-1, 1) + 1) / 2 

istart ( igrid, 1) = istart ( igrid-1, 1) + 

imax ( igrid-1, 1) *jmax ( igrid-1 , 1 ) * 
kmax ( igrid-1, 1 ) 


100 continue 


imax (1,2) = i 12 
jmax(l, 2) = jl2 
kmax(l,2) = kl2 

istart (1, 2) = istart (ngrid, 1) + imax(ngrid, 1) *jmax (ngrid, 1) * 

kmax (ngrid, 1) 


do 110 igrid=2, ngrid 

imax (igrid, 2) = (imax (igrid-1, 2) + 1) / 2 

jmaxdgrid, 2) = ( jmaxdgrid- 1,2) + 1) / 2 

kmax (igrid, 2) = (kmax (igrid- 1,2) + 1) / 2 

istart (igrid, 2) = istart (igrid-1, 2) + 

imax (igrid- 1,2) * jmaxdgrid- 1, 2) * 

kmax (igrid- 1,2) 


110 continue 


do 200 igrid=l , ngrid 

do 200 iblock=l,nblock 


62 



call metrics (x ( istart ( igrid, iblock) ) ,y ( istart ( igrid, iblock) ) , 
z ( istart ( igrid, iblock) ),... 
imax ( igrid, iblock) , jmax ( igrid, iblock) , 
kmax ( igrid, iblock) , . . ■ ) 

200 continue 


subroutine metrics (x,y , z , . . . , imax, jmax, kmax, . . . 
dimension x ( imax, jmax, kmax) ,y ( imax, jmax, kmax) , 
z ( imax , j max , kmax ) 


It is important to notice that in this example no changes were made to the subroutine 
metrics. Only changes to the main program were required and these involved the pointers. 
The important working subroutines in the computer code are the same for the multi- 
block program as for the multigrid computer code. New communication routines must 
be written to pass data between the blocks, but these are a separate issue and will be 
discussed in the Boundary Conditions chapter. 

Once the pointers are set up to accommodate a multi-block grid, then the DO loop 
structure must be modified to include loops to visit all of the blocks. The example for a 
W-Cycle has been modified to include multi-block flexibility and is shown by the bold 
characters in the following example: 


igrid=if ine ( 1 ) 

do 900 iblock=l , nblock 

call iterate { . . . , igrid, iblock, . . .) 
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do 5000 icycle=l , ncycle 
do 3000 ileg=l , nleg 


do 1000 igrid=ifine (ileg) +1, icoarse (ileg) , 1 

do 1000 iblock=l , riblock 


call restrict ( . . . , igrid-1 , igrid, iblock , . . . . ) 
call iterate ( . . . , igrid, iblock , • • - ) 


1000 continue 


do 2000 igrid=icoarse (ileg) -1, ifine (ileg+1) , -1 

do 2000 iblock=l , nblock 


call prolong ( . . . , igrid+1 , igrid, iblock , . . . ) 
call iterate ( . . . , igrid, iblock , . . .) 


2000 continue 
3000 continue 
5000 continue 
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At this point the advantages of grid topology independent subroutines has become 
apparent. The same subroutine can be used for all levels of multigrid for each block 
with no special logic required in the subroutine. Having only one subroutine to handle 
a specified set of tasks, regardless of the multigrid level on a block, also reduces the 
chances of coding errors, and the overall size of the source code. Also, a grid topology 
independent multi-block formulation allows the same executable to calculate the flows 
for a variety of geometric configurations. 

4.3 Multigrid Multi-Block Arrangements 

As stated in the multi-block structure section, there are two strategies that can be 
used in programming the multi-block — multigrid interaction. 

• Multigrid Inside of Multi-Block 

• Multi-Block Inside of Multigrid 

The first strategy is to have all multigrid levels run on a particular block before continuing 
on to the next block of the computational domain. This is referred to as the multigrid 
inside of multi-block method. This method lends itself to allowing an efficient use of 
memory space by computing one block at a time, writing that block out, reading in 
another block and operating on it before continuing to the next block. This is a very 
appealing way to handle computational problems that have large memory requirements. 
The disadvantage of this approach is that there is no communication between the coarse 
grid levels of the different blocks. At best this would only reduce the rate of convergence, 
because the lower frequencies would not be properly damped, if damped at all. At 
worst, the program will diverge, which becomes the case more times than not [68, 69]. 
Therefore, this approach was deemed unacceptable, because of its lack of robustness. 
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The second strategy, multi-block inside of multigrid, was chosen for this work. This 
approach operates with all of the blocks at the same multigrid level. Once the calculations 
for a particular multigrid level are completed, all of the blocks are processed to the next 
multigrid level, where operations are continued. A schematic of this approach is shown in 
Fig. 4.3.2. A problem with this approach is that if the memory requirements dictate that 
only one block be in use at a time, the input/output operation count becomes exceedingly 
high. The positive side of this approach is that this strategy mimics the same rate of 
exchange of information between cells as a single-block calculation. Thus, exactly the 
same convergence can be produced, and the original numerical efficiency of the single- 
block computer program is maintained, except for the time lost during the exchange of 
information across the block interfaces. The last statement is true only if the interface 
boundary conditions are of C 1 continuity, and there are no convergence acceleration 
techniques employed that are dependent on the placement of boundary locations. 


MULTIGRID /MULTIBLOCK STRATEGY 



Figure 4.3.2 Multigrid Multi-block Interaction Schematic. 
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4.3.1 Time Integration Strategy with Multigrid 


The calculation of a multigrid, four-grid, V-Cycle on a multi-block configuration, us- 
ing an n-stage Runge-Kutta time integration is outlined below. As previously mentioned, 
the type of multigrid multi-block interaction that was employed in the present work was 
multi-block inside of multigrid. The governing equation is 

N [Au‘] = -At[L'{u') - f] = RHS' (4.3.1) 

where N is the n-stage modified Runge-Kutta time integration scheme used to smooth out 
the high frequency errors, and i represents h, 2 h, 4 h, 8 h, etc. to indicate the different 
grid spacings. The numerical derivatives of the flux-vectors are symbolized by L'(u ), 
where l) is the operator providing the Euler or thin-layer Navier-Stokes flux derivatives, 
and u * represents the dependent variables. For the finest grid level, g h , f h is zero, because 
there is no multigrid forcing function for the finest grid level, only for the coarser grid 
levels. Starting the V-Cycle on the finest grid, g h , iterations are performed to smooth out 
the high frequency errors. One iteration involves computing all stages of the modified 
Runge-Kutta time integration method. The flux derivatives, which are called residuals, 
are computed for all blocks before continuing on to the next stage of the modified 
Runge-Kutta method. After each stage the dependent variables are updated, as well as 
the boundary and interface conditions. Also, if there is to be any residual smoothing, it is 
incorporated before the dependent variables are updated. After completing the necessary 
number of iterations to remove the high frequency errors, the dependent variables and 
the residuals are restricted to the next grid level, g 2h . Equation (4.1.15) provides the 
relation used to obtain f 2h . In this relation are two types of restriction processes. One 
is the restriction of the dependent variables ll h {u h ), and the other is the restriction of 
the residuals [ L h u h ] . The initial values of u 2h are obtained on the coarse grid by the 
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following volume weighted averaging, 


u 


2 h 



£ v<4 «f 

1=1 

£ v-oif 

«=1 


(4.3.2) 


The approximation sign is used because the error term, V 2h , shown in equation (4.1.16), 
is unknown and unaccountable (This initial value of u 2h will be used later as part of 
the correction for u h in the prolongation process from g 2h to g h ). The restriction of the 
residuals is performed by a simple summation of the residuals from the eight fine grid 
cells that compose a coarse grid cell, as shown in Fig. (B.1.4). The dependent variable 
residual restriction is given by 


8 


a“(«“) = /2‘[l‘.i‘] =E[ i * u 

1=1 


(4.3.3) 


With the initial u 2h values, the boundary and interface conditions are computed. The 
iterative process is performed again, eliminating the high frequency errors so that the 2 h 
grid data can be successfully restricted to the 4 h grid. The same equations that provided 
the relations for the transformation from g h to g 2h can be used to restrict data from g 
to g Ah . After the high frequency errors are eliminated on this grid level, the data can 
be restricted to the g u level. Finishing the necessary number of iterations on g 8h , the 
prolongation operation is then performed from g u to g Ah . The prolongation process 
provides a correction to the already existing u Ah values. This is done by calculating the 
difference between which was computed by the restriction process, and the u that 
has been updated by the iteration process. This is accomplished by the relation given in 
equation (4.1.17), which is rewritten here as follows, 

(““)«„ - (““L + «[«“ - *(““) J (4M) 

Then with the corrected u Ah more iterations are computed and then the correction data is 
prolonged to the next finer grid, g 2h . After iterating on this grid level the correction data 
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is finally prolonged back to the finest grid, g h . This completes one multigrid V-Cycle. 
The computation of a W-Cycle is similar. 


69 



Chapter 5 BOUNDARY CONDITIONS 


In general most CFD computer codes have what can be called standard inflow/outflow 
far field boundaries, inviscid wall, viscous wall, and symmetry plane boundary conditions. 
With a multi-block code there is an added interface boundary condition, which is used 
to exchange information between adjoining blocks. Each block is defined as a six sided 
volume, where each side is considered to be a face. For the present work, each face 
has two layers of ghost cell layers which contain boundary information. The ghost cells 
are used when the fluxes are being evaluated. Having two ghost cells allows higher 
order operations to be performed at block interfaces without any loss of accuracy or 
modifications. For solid wall boundaries only the ghost cell closest to the interior is 
used. It enforces no flow through the wall, and either parallel flow along an inviscid wall 
or no slip on a viscous wall. The only contribution a solid wall provides to the flux is 
the pressure term in the momentum equations. For the symmetry plane, inflow/outflow 
far field boundaries, and interfaces, the fluxes are computed at the boundary in the same 
manner as for an interior point, because both layers of ghost cells on the block face 

are used. 

Boundary conditions normally drive the flow field. Ghost cells are used as an aid 
in enforcing boundary conditions accurately. These cells can be used in different ways. 
One way is to treat them as storage space to maintain dependent variable values on a 
wall. They are also used to maintain the far field boundary information, which controls 
the type of flow field the geometrical configuration will encounter. 
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The present control volume approach uses ghost cells as an extension of the interior 
cells of the computational domain. This allows the ghost cells to have assigned volume 
and area vectors. For external far fields, the ghost cell values are determined based on 
the prescribed flow conditions and interior cell values of the computational domain. The 
intent is to have ghost cell values that represent what the correct flow values would be 
in those specific x, y, and z locations for an identical physical flow field. 

In treating the wall boundaries, the ghost cells again have volumes and area vectors, 
but rather than store flow information on the wall, they maintain values that when 
combined with the first interior cell, normal to the wall, provide the correct wall flow 
quantities. This approach is useful in evaluating the fluxes, because it allows the gathering 
of cell information for the first interior cell, normal to the wall, to be conducted with the 
same approach as any other interior cell in the domain. 


5.1 Far Field Inflow/Outflow Boundaries 

For the far field inflow/outflow boundaries the Mach number normal to the cell face 
is computed to determine if the flow is subsonic or supersonic at the inflow/outflow 
boundary. If the flow is supersonic then the direction of the flow is determined and a 
direct transfer of the reference conditions to the ghost cells is used for the inflow case. 
For outflow, a linear extrapolation of the interior values to the ghost cells is used. If the 
flow is subsonic then the one-dimensional characteristic equations are used to determine 
the ghost cell values. The two Riemann invariants, R + , and R~, are computed, and their 
average taken to give the velocity normal to the cell face, q n , which is used to indicate 
whether the case is an inflow or outflow boundary condition. From the one-dimensional 
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characteristic equations, the following formulations are obtained, 

Qref = u refh + V re fl y + W Te fl 2 

Qint = Umtlx + V int l y + W xnt l z 

, 2 
I\ — Qint H -(lint 

7 - 1 
2 

R — Qref ~ j a ref 


( 5 . 1 . 1 ) 


Qn = 7 >(^ + + R ) 


/ — — where / = £, 77, & an< ^ w — Xi Vi & z 

y/P, + ? + 'I 

with a being the speed of sound, and subscripts re/ and inf indicating reference and 
interior, respectively. If the normal velocity, q n , is negative, then the boundary condition 
is subsonic inflow. The ghost-cell, Cartesian velocities are computed based on the 
reference Cartesian velocities, u re/ , »„/, and w Tef , and the difference between the 
average Riemann normal velocity, q n , and the reference normal velocity, q ref as follows; 


5 * = Pref / P re f 

U g = U Te f + (q n ~ qref ) lx 

V g = V Te f + (q n ~ q r ef)ly 
W g = W re f + (q n - qref)h 


( 5 . 1 . 2 ) 


where s* is an isentropically derived entropy value, and g denotes the ghost cell values. 
This formulation stipulates that the entropy will not change across the far field boundary. 


For the subsonic outflow case; 


s * — Pint / P* n< 

Ug = Ui n t + (<?n — q\nt)lx 
V g = V ln t + (q n ~ qint)ly 

Wg : W tn i + {q n — qtnt)lz 
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The density and total energy per unit volume are then computed as follows, 

«= j(7-l K« + -*') 

__ / a 2 s*\ ~ 

P9 ~ \~T ) (5.1.4) 

_ a 2 pg 

p 9 7 

E g = : ^ + \p g {u 2 g + v 2 g +w 2 g ) 

5.2 Symmetry Plane and Solid Wall Conditions 

For the symmetry plane boundary condition, the velocity vectors of the two interior 
cells adjacent to the symmetry plane are reflected into the ghost cells, and the densities 
and total energies per unit volume are transferred directly. The inviscid wall conditions 
reflect the velocity vector of the first interior cell across the wall surface, maintaining 
the tangential velocity, just as with the symmetry condition. The pressure and density 
are obtained from the interior cells by either a direct transfer or a linear extrapolation. 
The total energy per unit volume is then computed using the ghost cell values. For 
the viscous wall the same procedures are used with the exception of having a no slip 
condition on the wall. 

5.3 Block Interface Conditions 

The interface conditions require that the blocks meet with C 1 continuity, meaning 
that the grid lines are continuous and that continuous grid metrics are maintained across 
block interfaces. This approach allows ghost cells of one block to receive information 
by direct transfer from the corresponding cells in the adjoining block; therefore the 
process is a one to one transfer of data, with no averaging or approximation. Having 
the two layers of ghost cells on each block face allows operations at the interface to be 
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performed in the same manner as any operation would be performed in the interior. There 
is no issue concerning the conservation of mass or a flux balance across the interface 
due to the continuity requirement across the interface. Plus, this approach provides 
the same order of accuracy at the interface cells that is present in the interior cells. 
All other interface conditions cannot provide both higher order accuracy and maintain 
a complete conservation of the fluxes across an interface. One advantage of this type 
of interface condition is that it allows the computer code to obtain exactly the same 
convergence for a case that could be run as a single block. This was very helpful as a 
debugging tool during the initial development of the computer program. The computer 
code was developed to allow any face of one block to interface with any face of another 
block. This required that the indices of one block be allowed to adjoin different index 
families of the adjacent block. Also the interface routine permits the indices to increase 
in the same direction or in opposite directions across a block interface. All that is 
required is that both blocks adhere to the right hand rule. More sophisticated interface 
conditions can be implemented by changing the interface routine to one that best suits 
a particular configuration’s requirements. The rest of the computer code will remain the 
same, because at no other time is it necessary for the blocks to interact with each other. 
That is why the two layers of ghost cells on each block face are used. Once these cells 
have been updated the blocks are not required to communicate with each other until the 
next iteration. Each block is self-contained, allowing it to go through the flux evaluator 
and the prolongation and restriction multigrid routines without requiring any additional 
cell information from the other blocks. 

Each face of a block can be divided into multiple segments or patches, and each 
patch can have a different boundary condition. The only constraint on these patches is 
that they need to be of multigridable indices, so that the physical x, y, and 2 locations of 
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the patch will not change when the solver is on a different multigrid level. This is very 
important, especially if one of the multiple patches is interfacing with another block. 


5.4 Definition of Multigridable Index 

A multigridable index is an integer value based on 2 n +l. An example of a multi- 
gridable index is the number 129, which is 2 7 +l. It is multigridable because if this 
were the number of points in a particular coordinate direction, and every other point was 
eliminated, the total number of points would become 65, which is 2 6 +l, and the 65 th 
point would be in the same x, y, and z location as the 129 th point of the finer grid. Also, 
if every other point were eliminated again, the total number of points would become 33, 
which is 2 5 +l. The 33 rd point would also be in the exact x, y, and z location as the 65 th 
and the 129 th points of the finer grids. This would not be the case if the initial number of 
points on the fine grid was 1 12. If every other point is eliminated, starting from the first 
point, the total number of points would be 61, and the 61 st point would not be in the same 
x, y, and z location as the 112 th point. This is very important when multiple boundaries 
are used, because when patches are involved, the starting and ending index needs to be 
multigridable; otherwise error could be introduced into the calculations due to moving 
boundary locations or the exchange of incorrect information across interfaces. When the 
restriction process to remove every other point is executed on a 65-point grid, every value 
less than 65 that is a multigridable index will maintain the same x, y, and z location. 

It should be noted however, that the convergence rate for a problem can change 
if a single-block domain is divide into multiple blocks when residual and/or corrector 
smoothing is being employed. The smoothers employed in the present computer code only 
operate on one block at a time, and if a domain has been divided into smaller segments, 
then the influence of boundaries at one end of the entire computational domain will take 
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longer to reach flow field data at the other end of the domain. For subsonic flows this 
could definitely slow the rate of convergence, whereas it may be beneficial to supersonic 
flows, based on the physical directions of information propagation. Therefore, dividing a 
single block into multiple blocks may provide better or worse convergence, depending on 
the physics of the flow when various domain boundary-dependent acceleration techniques 

are employed. 
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Chapter 6 CASES STUDIED 


The flow cases presented here were chosen to validate the computer code, and display 
its flexability in accommodating different types of geometrical configurations. Known 
flow cases were chosen to validate the incorporated flow solvers. This was done to 
insure that the employed acceleration techniques of residual and corrector smoothing, 
and multigrid acceleration were not introducing error into the final numerical results. 
The flexibility was examined by testing a variety of geometric flow configurations. 
First a comer flow case was computed using a single block, then using a multi-block 
configuration. Next the multiple inflow boundary condition was tested, followed by a 
multiple boundary condition including a block interfacing with itself. The final case 
tested required true multi-block capabilities, coupled with the ability to accommodate an 
interface between blocks with different mesh topologies. All of the flow configurations 
were computed with the same computer code, only changes to the input file were required; 
therefore proving the computer code’s flexibility. 


6.1 Inviscid Corner Flow 

The first case studied was an inviscid corner flow. The main reason this case was 
chosen was to demonstrate that the present multi-block method allowed the interior cells 
at the interfaces of the blocks to be treated as any other interior cell. This approach 
allowed the multi-block configuration to produce the exact same results as the single-block 
configuration. This was demonstrated by computing the flow field using a single-block 
grid configuration, and then dividing that single block into eight blocks and computing 
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the flow field again. By doing so, this case showed what is required of the interface 
routine for similar flow configurations. Finally, the three-dimensional flow effects of this 
case required that there be no coordinate direction biasing of the flow solution by the 
computer code. 

A compression comer was generated by placing a compression ramp on the bottom 
and back walls of a rectangular duct. A schematic of the comer generated from the 
connection of the back and bottom walls is shown in Fig. 6.1.1. Indicated in the figure 
are Pb, which was taken as a pressure reference point, and Y c , which indicates the y 
location at the intersection of the compression ramps. These values were used m results 
comparisons. The supersonic inlet flow was M tnkt = 3.0, and the ramp angles were 
a = 9.5°. Thus, the grid and therefore the flow was symmetric about the intersection 
of the back and bottom walls, forming the compression comer. For supersonic corner 



Figure 6.1.1 Schematic of Compression Comer Duct. 
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flow, three shock structures are produced. Two of the shock structures are wedge flow 
shocks, and their flow properties can be verified using two-dimensional analysis based 
on the Mach number normal to the leading edge of the wedge. The third structure is 
produced when the two wedge shocks coalesce to form a three-dimensional flow region 
shaped like a cone, with its apex located at the intersection of the leading edges of the 
compression ramps. The base of the cone is coincident with the exit plane of the duct. 
These characteristics can be seen in Fig. 6.1.2, where Mach line contours are projected 
on the back and bottom walls, and on the exit plane of a 49x49x49 channel or duct gnd. 
This case was performed using FMG, no flux limiter, and modified Runge-Kutta time 
integration. The position of the wedge shocks is shown on the back and bottom walls 
by the region of highly concentrated Mach lines perpendicular to the inflow direction. 
The edges of the cone shaped shock surface can also be seen on these two walls. Four 
flow regions are present on the exit plane. In the upper right comer is free-stream flow, 
which is one-dimensional. From the middle of the plane to the lower left comer, the 
flow is three-dimensional, and the bottom of the cone surface, a partial disc, can be seen. 
The wedge shock planes, which are two-dimensional flows, can be seen in the upper left 
corner and the lower right corner of the exit plane. Note that since the geometry of the 
channel is symmetric about the compression corner (the one joining the back and bottom 
walls to the exit plane of the channel), the flow field should be and is symmetric about 
this corner. Corner flow shock structures are identified by having triple points, where the 
three-dimensional flow region meets the wedge shock flow and the one-dimensional free- 
stream flow region. The results shown here have two triple points. Shown in Fig. 6.1.3 
is a schematic of how a plane parallel with the exit plane of the channel would appear. 
It shows the flow structure that results from this type of comer flow (figure based on 
information obtained from Ref. [98]). The present flow results are in good agreement 
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Figure 6.1.2 Mach Line Contours, M inlei = 3.0 and a - 9.5°. 


WALL 

SHOCK 



with the numerical results of Marconi [98] and Kutler [99]. Shown in Fig. 6.1.4 is 
a comparison of the present work with the shock fitting numerical results of Marconi 
[98], the finite difference results of Kutler [99], and the experimental results of West and 
Korkegi [100]. It shows the relative pressure distribution on the surface of the back wall 
along a coordinate line that is parallel to the exit plane of the channel (Note that since the 
geometry of the channel is symmetric about the compression comer, either the back or 
bottom wall could have been used). The results of Marconi and Kutler were plotted from 
data presented in Ref. [98], where they solved for the comer flow in a two-dimensional 
plane, similar to what is shown in Fig. 6.1.3. The range of the pressure distribution for 
the present work is very similar to Kutler’ s predictions, but the shock is smeared in the 
present study because its structure is skewed relative to the computational grid, and the 
grid spacing in the streamwise direction is rather coarse. 

To better capture the shock, the same test case was performed on a 65x65x65 grid. A 
comparison between the 49x49x49, the previously mentioned grid and a 33x33x33 grid 
is shown in Fig. 6.1.5, where it can be seen that the shock is less smeared on the finer 
grids. The packing of points in the streamwise direction has a very significant influence 
on this type of results comparison. Convergence histories are shown in Fig.6.1.6 for this 
comer flow case using different extrapolation techniques. As expected, the first-order 
extrapolation converged the quickest, followed by Fromm’s method and k= 1/3. The 
residual for the pure second-order upwind method became hung after converging nine 

orders. 
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Figure 6.1.4 Comparison of Relative Pressure Distributions on Wall of Corner Flow. 
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Figure 6.1.5 Grid Refinement Study on Comer Flow. 
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Figure 6.1.6 Convergence Histories for the Comer 
Flow Using Different Extrapolation Techniques. 
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6.1.1 Multiple Interface Requirements 


To test the multi-block capabilities, the 65x65x65 comer flow grid was divided into 
eight blocks by bisecting each coordinate direction. A schematic of the block domains 
is shown in Fig. 6. 1.1.7, where the shaded surfaces indicate the interfaces between the 
blocks. For computing the fluxes, the exchange of information across an interface needs 



Figure 6. 1.1.7 Eight-Block Configuration for Compression Corner. 

to include only the ranges of the interior cells on that interface. However, when using 
multigrid techniques, consideration must be given to the restriction and prolongation 
operations. Whether or not there are interfaces has no influence on the present restriction 
process, but it can have an influence on the prolongation process. By dividing the 


85 



domain into blocks there are block boundaries where there would have normally been 
standard interior cells. This can affect the values produced by the trilinear interpolation 
prolongation process. The prolongation process is explained in Appendix B. The standard 
procedure for the present prolongation process at a boundary is to set the corrections to 
zero. That is because the boundary ghost cells are not solved directly; they are updated 
based on the flow values of the interior cells. This indicates that the correction at a block 
interface should be treated differently than that of a solid surface. Interface corrections 
should allow the prolongation process to compute the same results for a divided region 
as it would for an undivided region. To properly execute the prolongation process at an 
interface requires computing the correction for the ghost cell values and using them in 
the same manner as the correction values of the interior cells. The problem arises at the 
edges of the interface. It is these locations that require information from the ghost cells 
normal to the block interface, and also from the cells diagonal to it. This can be seen 
in Fig. 6. 1.1.8, where a plane from four interfacing blocks is shown. The thick lines 
represent the interior regions of the blocks, and the thin lines represent the ghost cells. 
If the four blocks were joined as one, the prolongation for the upper left corner of block 
1 would be influenced by information in the other three blocks; therefore the range of 
information transferred should fill the ghost cells of block 1 with the first two bottom 
rows of “x’s” from block 2, the four lower right corner “o’s” from block 3 and the first 
two right-hand-side columns of “z’s” from block 4. Transferring these values then allows 
the prolongation for the upper left corner of block 1 to produce the same values as would 
be produced if the four blocks were combined as one. It is important to note that block 1 
interfaces with only blocks 2 and 4; therefore, the range of information gathered across 
an interface needs to include more than just the values from the interior cells; it needs to 
obtain two extra cells from each end of the interface. A range which includes two extra 
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Interior Cells 




cells at each end of the interface provides the “o”s from block 3 to block 1, and allows the 
ghost cells of each block to have all of the information necessary to correctly complete 
the prolongation process. Since information is exchanged between at most two blocks at 
a time, two passes through the interface routine are required so that all of the interface 
ghost cells (beyond the interior index range of the interface) will be updated with the 
current interior cell values from the different blocks. This is necessary since the interface 
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ghost cells are updated in a sequential manner, generally starting with block 1 interfaces, 
then proceeding to block 2 and so on. As one can see from Fig. 6.1.1.8, the interface 
ghost cell values required by block 1 from block 3 will lag because block 1 gets its final 
block 3 information from block 4. Therefore, until block 4 has had its interfaces updated 
with current block 3 information, the information block 4 passes to block 1, on the first 
pass through the interface routine, will be from the previous time step. This problem is 
eliminated by executing the interface routine twice. Also, this method requires that block 
1 only know the blocks that are normal to its interfaces, and not what blocks are diagonal 
to its comers, which in this test case is block 3. This simplifies the input information 
about the boundary conditions. Also, this interface approach satisfies the requirements 
of a block having three connected surfaces as interfaces, which occurs in this eight-block 
corner flow configuration. Notice that every block in this configuration interfaces with 
three other blocks, as can be ascertained from Fig. 6.1. 1.7. 

In Fig. 6. 1.1. 9 and Fig. 6.1.1.10 the Mach line contours from the single-block and 
eight-block calculations are shown, respectively, on the back and bottom walls, and the 
exit plane of the channel. For the eight-block configuration, the interfaces between the 
blocks are shown as thick solid lines on the back and bottom walls, and exit plane. As 
can be seen, the Mach lines pass through the interfaces without deviation, proving that 
the interfaces are not introducing any error. Similar results are displayed in Figs. 6. 1.1.1 1 
and 6.1.1.12, where the pressure contours are shown for both the single- and eight-block 
configurations, respectively. Again, there is no difference in the contours between the 
single-block and the eight-block calculations. As shown in Fig. 6.1.1.13, the convergence 
histories for the single-block and the eight-block calculations are the same. 
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Figure 6. 1.1.9 Mach Line Contours for the 
Single-Block Calculation, M, n i et = 3.0 and a = 9.5°. 



Figure 6.1.1.10 Mach Line Contours for the 
Eight-Block Calculation, M in i e t = 3.0 and a = 9.5°. 
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Figure 6.1.1.11 Pressure Line Contours for the 
Single-Block Calculation, Mi n i et = 3.0 and a = 9.5°. 



Figure 6.1.1.12 Pressure Line Contours for the 
Eight-Block Calculation, M iniet = 3.0 and a = 9.5°. 
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Figure 6.1.1.13 Comparison of Convergence Histories Between the 
Single-Block and Eight-Block Calculations for the Corner Flow. 
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6.2 Pseudo Two-Dimensional Jet Exhaust Plume 

The second case studied was an inviscid pseudo two-dimensional jet exhaust plume. 
This case tested the multiple boundary input format of the computer code as well as 
provided a basis for comparing the two different upwind solvers, in terms of their ability 
to predict both slip lines and shocks. It also served as a forum for the examination of 
the different extrapolation techniques of pure second order upwind, Fromm s scheme, 
and k = 1/3. This case is considered a pseudo two-dimensional flow because the cross 
flow is negligible. 

The inflow surface of the computational domain was divided into two sections. One 
section had as its inflow conditions the exhaust from a jet of height “h , where M }e t = 1.5. 
The second section had a free-stream inflow of Moo = 2.5. The static pressure and 
static temperature ratios between the two flows were p je t/Poo - 3.5 and T je tjT 0 © 3.0, 

respecitvely. A schematic of this flow field is provided in Fig. 6.2.14. This figure 
identifies the slip line, expansion fan, and shock that is present in this type of flow. 
The “wall” in this figure was actually treated as a symmetry plane. The grid generated 
for the pseudo two-dimensional plume flow is shown in Fig. 6.2.15. This case was 
first tested incorporating Roe’s flux-differencing using: (1) pure second order upwind, 
(2) Fromm’s scheme, and (3) k = 1/3 extrapolation of the primitive variables, without 
the use of a flux limiter. These results were compared with those from a validated shock 
fitting code developed by Salas [101]. This case was evaluated using FMG, variable 
coefficient residual smoothing and modified Runge-Kutta time integration. At the first 
station, x/h « 1.0, shown in Fig. 6.2.16a, the second order extrapolation produced a 
larger undershoot at the slip line than the other higher order methods. It is comparable 
everywhere else, produces a smaller overshoot at the shock, and has no oscillations. The 
k = 1/3 continually produces the largest overshoot at the shock and is accompanied by 
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a small oscillation (Figs. 6.2.16d - 6.2.16e), up until station x/h « 5.0 (Fig. 6.2.16e). 
At this station and station x/h = 10., shown in Fig. 6.2. 17f, the overshoots at the shock 
are greatly diminished for all extrapolation methods. The two vertical columns along 
the right side of the figures indicate the packing density of points in the y/h direction. 
The shock fitting computer code generates its own grid, and its packing density is shown 
as the first column starting on the left. The Roe scheme employed a different grid, 
which had a vertical packing density shown as the right column. The columns are shown 
in the order of the flow solver legend, with the far left column representing the grid 
spacing for the last flow solver in the legend. Based on these results, and the desire 
not to employ a limiter (because of the convergence difficulties they produce), the pure 
second order upwind extrapolation method was chosen for the next case. Convergence 
history comparisons between the different extrapolation techniques for the pseudo two- 
dimensional jet exhaust plume are shown in Fig. 6.2.17. 


Y 



Figure 6.2.14 Schematic of Pseudo Two-Dimensional Jet Exhaust Plume. 
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Figure 6.2.15 Grid Used for the Present Pseudo 
Two-Dimensional Jet Exhaust Plume Calculations. 


The next comparison for this case was between Roe’s flux-difference splitting and 
van Leer’s flux-vector splitting. Figures 6.2.18 and 6.2.19 show the Mach line contours 
for Roe’s flux-difference splitting method, and van Leer’s flux-vector splitting method, 
respectively. In Figs. 6.2.20 and 6.2.21 are the pressure contours for Roe’s scheme and 
van Leer’s scheme, respectively. Both methods employed the same Runge-Kutta time 
integration scheme, using pure second order upwind extrapolation, without flux limiters, 
and they employed variable-coefficient residual-smoothing. Conservative variables were 
extrapolated for van Leer’s scheme, because the use of primitive variables is known to 
cause oscillations. Primitive variables were extrapolated for Roe’s scheme. The Mach 
line contours appear to provide the same solution, with van Leer’s method spreading the 
slip line a little wider than Roe’s method. This is to be expected, because flux-vector 
splitting does not have a mechanism for resolving the slip line, which flux-difference 
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splitting does because it is based on an approximate Riemann solver. Overall the 
differences in the two sets of contours are negligible. Further comparisons of these 
results are shown in Figs. 6.2.22a - 6.2.22f, with the validated shock fitting code of 
Salas [101]. The first figure (Fig. 6.2.22a) is for station x/h « 1.0, the only disparity 
between the flux-vector splitting and the flux-difference splitting methods comes just 
after the constant Mach number region and just before the slip line. In that region, the 
flux-difference splitting method has a larger undershoot Both upwind methods have the 
same magnitude of undershoot and overshoot at die shock, and neither had an overshoot 
at the slip line. The shock fitting method tends to have oscillations in the transition 
region between the expansion fan and the constant Mach number region, but provides 
crisp results throughout the remainder of the flow for this x/h location. For x/h ~ 2.0, 
shown in Fig. 6.2.22b, one can again see that the only disparity between the van Leer 
and Roe methods is the undershoot by Roe’s scheme at the slip line. Starting at this 
x/h location, the overshoot at the shock has essentially disappeared. This is due to the 
shock becoming aligned with the grid. The next location, x/h ss 2.5 (Fig. 6.2.22c), 
shows that the undershoot from Roe’s scheme at the slip line has diminished. At the 
x/h w 3.0 location, shown in Fig. 6.2.22d, the Roe and van Leer methods are providing 
the same magnitude of undershoot at the slip line, but at x/h k 5.0 (Fig. 6.2.22e) and 
x/h = 10.0 (Fig. 6.2.22f) the two upwind methods are providing essentially the same 
error in the overshoots and undershoots; however, Roe’s scheme again gives a slighdy 
larger undershoot at the slip line. 

Overall the agreement between the two upwind methods was very good and their 
results compare well with those of the shock fitting computer code of Salas [101]. Both 
Roe’s and van Leer’s upwind methods predict overshoots and undershoots at slip lines 
and shocks on a grid that is not properly aligned with the flow physics. 
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Figure 6.2.16 Comparisons Between Different Extrapolations of Roe’s Scheme and a 
Shock Fitting Code for a Pseudo Two-Dimensional Exhaust Plume. (Continued . . . 
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Figure 6.2.16 Comparisons Between Different Extrapolations of Roe s Scheme 
and a Shock Fitting Code for a Pseudo Two-Dimensional Exhaust Plume. 
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Figure 6.2.17 Convergence History Comparisons Between Different 
Extrapolation Techniques for the Pseudo Two-Dimensional Jet Exhaust Plume 




Figure 6.2.18 Mach Line Contours for Roe’s Scheme, 

Moo = 2.5, Mjet = 1.5, Pjet/Poo = 3.5, and T je t/T x = 3.0. 



Figure 6.2.19 Mach Line Contours for van Leer’s Scheme, 
Moo = 2-5, M jet = 1-5, P je t/Poo = 3.5, and = 3.0. 
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6.3 Laminar and Turbulent Flows Over a Flat Plate 


The third test case was laminar and turbulent flow over a flat plate. The purpose 
of this case was to demonstrate that the viscous terms had been implemented correctly 
and that the Baldwin-Lomax algebraic turbulence model had been employed correctly 
for attached turbulent flows. 

6.3.1 Laminar Flow 

For laminar flow over a flat plate, the Mach and Reynolds numbers were taken to 
be Moo = 0.5 and R = l,000/(unit length), respectively. The normalized height of 
the first cell normal to the plate was lxlO -4 units, with one unit being the reference 
length of the plate. Good agreement was obtained between the present results and the 
Blasius solution. Comparisons with the skin friction coefficient and velocity profile 
are shown in Fig. 6.3.2.23. The grid was 65x65x5, with 64 cells in the streamwise 
direction, 64 cells normal to the plate, and 4 cells in the span-wise direction, to allow 
for three levels of multigrid. Computations were performed incorporating FMG, variable 
coefficient residual smoothing, and modified Runge-Kutta time integration. There was a 
minimum of 34 grid cells in the fully developed boundary layer (where fully developed 
is defined as having self similar velocity profiles). Convergence histories for the laminar 
flat plate flow comparing different acceleration techniques are shown in Fig. 6.3.2.24. 

As one can see, multigrid acceleration coupled with residual smoothing provided the best 

convergence rate. 

6.3.2 Ttirbulent Flow 

The turbulent flat plate flow conditions were M = 0-5, R = 1, 000, 000/(umt length), 
and the normalized height of the first cell normal to the plate was lxlO" 5 units, with 
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one unit being the reference length of the plate. Good agreement was obtained for the 
skin friction, velocity profile, and law of the wall plot, as shown in Fig. 6.3.2.25. The 
grid was 65x95x5, with 64 cells in the streamwise direction, 94 cells normal to the plate, 
and 4 cells in the span-wise direction, to allow for three levels of multigrid. Compu- 
tations were performed incorporating FMG, variable coefficient residual smoothing, and 
modified Runge-Kutta time integration. A minimum of 42 grid cells were in the fully 
developed boundary layer (fully developed meaning self similar velocity profile). 




Figure 6.3.2.23 Laminar Flat Plate Comparisons with Analytical Calculations 


106 


Log (Residual) 



Figure 6.3.2.24 Convergence Histories for the Laminar Flat 
Plate Flow Comparing Different Acceleration Techniques. 
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Figure 6.3.2.25 Turbulent Flat Plate Comparisons with Analytical Calculations 



6.4 Turbulent Flow Over an ONERA M6 Wing 

The fourth case studied was turbulent flow over the ONERA M6 wing [102] with 
193x49x33 grid points in a C-0 mesh topology (schematic shown in Fig. 6.4.26). 
This case was chosen for two reasons. First, this case has true three-dimensional 
turbulent flow. Second, this flow configuration places special requirements on the block 
interface routine, which will be explained in the next section. The first test case was 
Moo = 0.699, a = 3.06°, and R = 11.7xl0 6 /(imit length). The wing was normalized 
to a semi-span of a unit length. This was a subcritical case shown to valid the computer 
code, which is in good agreement with the experimental data [102], as indicated by the 
C p plots in Figs. 6.27a - 6.27f, where t] is the dimension distance from the wing root. 




Figure 6.4.26 Schematic of C-O Mesh Topology for ONERA M6 Wing. 
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Figure 6.27 Comparison of Numerical Results 
with Experimental Data for ONERA M6 Wing. 
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The second test case has a lambda shock structure on the upper surface, and 
is much more demanding numerically, with Moo = 0-84, a = 3.06°, and 
R = 11.7xl0 6 /(unit length). Since the Reynolds number is the same as the previous 
case, the same grid was used. Note that a significantly higher Reynolds number would 
require smaller vertical spacing between the points in the boundary-layer regions. Com- 
parisons of C p plots with experimental data and other numerical results are shown in 
Figs. 6.4.28a - 6.4.28f. TLNS3D is a Jameson type, thin-layer, central difference, modi- 
fied Runge-Kutta computer code, developed by Vatsa and Wedan [86], which employed 
the Baldwin-Lomax algebraic turbulence model [72]. CFL3D is an implicit approximate 
factorization, thin-layer, upwind computer code using Roe’s flux-difference splitting, and 
it also utilized the Baldwin-Lomax algebraic turbulence model [72], This computer code 
was developed by Thomas, Krist, and Anderson [103]. The present results used Roe’s 
flux-difference splitting, with pure second order upwind extrapolation, which did not 
require a flux limiter. As can be seen in these figures, the present results are in good 
agreement with the experimental data [102], except at the second shock location at the 
rj = 0.80 station, shown in Fig. 6.4.28d. Here all of the numerical results disagreed 
with the experimental data. 

The convergence history for this case can be seen in Fig. 6.4.29. Although the 
residual, which is the L 2 norm of the density, only converges three and a half orders, 
the coefficients of lift and drag and the number of normalized supersonic points are 
converged in less than two hundred work units. 
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(a) 



(b) 

Figure 6.4.28 Comparison Between the Present Results and Other 
Numerical Results for the ONERA M6 Wing. (Continued . . . ) 
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(C) 



(d) 

Figure 6.4.28 Comparison Between the Present Results and Other 
Numerical Results for the ONERA M6 Wing. (Continued . . . ) 
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Figure 6.4.28 Comparison Between the Present Results 
and Other Numerical Results for the ONERA M6 Wing 
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Figure 6.4.29 Convergence History for ONERA M6 
Wing at A/oo = 0.84, a = 3.06°, R = 11.7 x 10 6 /imit. 
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6.5 Requirements for Non-Interface with Interface 
Multiple Boundary Conditions 

The block interface difficulty comes at the intersection of the trading edge of the 
wing and the wake region. When using one block to solve for the wing flow, the wake 
region is treated as an interface. In this case the “j = jmin” face of the block, which is 
a constant ( surface, has three boundary conditions. Starting at “1 = min”, which is the 
lower half of the exit plane, and continuing to the trailing edge of the wing, the jmin 
surface boundary condition is an interface across the wake cut. Around the wing, the 
boundary condition is a turbulent solid wall. Continuing from the upper wing surface 
trailing edge to “imax”, which is the upper half of the exit plane boundary, the “jmin” 
surface boundary condition is an interface across the wake cut to the first “jmin” interface 
boundary condition; therefore across the wake the “jmin” face of the block exchanges 
information with itself. The interface conditions used for the comer flow test case, where 
the grid was divided into eight blocks, will cause problems at the trailing edge of the 
wing. The process of gathering variable information from the extra ghost cells across the 
interface becomes fatal in this case because it replaces the solid wall boundary conditions 
with interior flow cell values. This occurs at the first two pairs of cells, at the trailing 
edge of the wing on both the upper and lower surfaces. This can be seen in Fig. 6.5.30, 
where the affected cells are indicated by the octagons and circles. The thick lined region 
represents the interior cells, and the thin lined region represents the ghost cells. The 
surface between the two regions is part solid wall, which is the wing indicated by the 
hash marks, and the dotted lines indicate the wake cut. If the interface conditions were 
not adjusted for this particular case, the values for the ghost cells on the lower side of the 
wing would be the “y’s” in the octagons rather than the “sL” values, which are the solid 
wall boundary conditions for the lower surface. The same problem exists for the circled 


118 


ghost cells on the upper side of the wing. Obviously these ghost cells need to be set to 
their proper values before they are used. This is a problem for this particular situation 



only because a non-interface and an interface boundary condition are both used on the 
same block face. This problem can be eliminated by enforcing the solid wall boundary 
condition again, after the interface routine has been executed. To identify this problem 
requires an evaluation of the types of multiple boundaries prescribed for a block face. 
Note that at the “kmax” face, which is at the wing tip and is where the C-0 mesh folds 
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or closes onto itself, the interface condition is used again to allow the “kmax” surface to 
exchange information with itself. The bottom half of the C mesh, which starts at “imin” 
and continues to the leading edge of the wing, exchanges information with the top half 
of the wing which starts at the leading edge and continues to “imax”. 

6.6 Afterbody with Internal Nozzle 

The final test case required a multi-block code to handle internal and external flow 
for an afterbody configuration, shown in Figs. 6.6.31 and 6.6.32. This afterbody 
was examined in a wind tunnel at various free-stream Mach numbers by Putnam and 
Mercer [104], and Compton, Thomas, Abeyounis, and Mason [105]. For certain cases 
C p data were obtained from the center of the top and side of the boat-tail. The external 
geometry was relatively easy to grid, except for the boat-tail, which required the use ot a 
ninth order super elliptic equation, with various offsets to produce the correct curvature 
for the boat-tail comer edges. The formulations and geometry of this afterbody are 
described in reference [104], In the present work a polar grid was used to generate the 
outer surface geometry of the afterbody. This can be seen in Fig. 6.6.33. A closer 
view of the afterbody is shown in Fig. 6.6.34, displaying some its grid structure, where 
only everyother point is shown on the afterbody for visual clarity. Using a polar grid 
allowed for the direction normal to the external surface to always be in the r? direction. 
This in turn provided the necessary direction for the length scale that was needed for the 
algebraic turbulence model. The more demanding part of the grid generation process was 
the interior nozzle. Air was supplied subsonically at a specified temperature and pressure 
into a circular cross sectional area (section FS 40.95 Fig. 6.6.31). The flow passed through 
a baffle plate and then continued toward a settling chamber, with the cross sectional area 
changing smoothly from circular to rectangular. The width of the internal geometry 
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remains constant, as the upper and lower surfaces converge to produce a sonic throat 
region. The nozzle re-expands producing a supersonic exit flow of M jet exit = l- 6 - Th e 
nozzle exit is rectangular, which is the type of exit a rectangular thrust vectoring/thrust 

reversing nozzle would have. 

The grid for the internal nozzle was very demanding. First, the initial cross section 
was circular, which was best suited by a polar grid; however, once the cross sectional area 
changed to rectangular, a polar grid was no longer appropriate. In fitting a polar grid to a 
rectangular cross section, it becomes extremely difficult to force the radial lines, extending 
from the polar axis, to be normal to the sides of the rectangular cross section. Also, many 
computer codes have difficulty accommodating singular grid lines, such as the polar axis 
in this case. The problem generally stems from the fact that at singular grid lines, one 
face of the control volume cell has a zero area. This problem is accentuated when the 
grid lines are packed in the circumferential direction, creating highly skewed cells. One 
way to alleviate both problems is to use an H-H grid topology for the internal geometry. 
This eliminates the singular line problem and the grid lines are naturally normal to the 
solid surfaces. To capture the turbulence effects at the walls requires a dense packing 
of points, which with the H-H topology results in an extremely large number of points 
in the comers of the rectangular cross section. Also, at the exit of the nozzle, which is 
at the end of the afterbody, the grid lines from the internal grid are required to match 
the external grid lines with C 1 continuity. This is an obvious problem. For the external 
grid to have such a large cluster of points at the comer edges of the afterbody would 
have required a tremendous number of grid lines, which is infeasible with the current 
memory restrictions on today’s super computers. Furthermore, such a large number of 
points would require a large amount of CPU time to achieve converged flow solutions. 
The best approach for griding the interior nozzle is to use a polar grid at the solid wall 
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surfaces, pack the necessary number of points for the turbulent boundary layer, and then 
have the polar grid interface with an H-H grid, which then fills in the remainder of the 
interior of the grid and eliminates the singular line. Three a— locations wre chosen to 
show the interfacing of the H-H topology and the polar topology in Fig. 6.6.35. In this 
figure, every third point is shown for the polar topology for visual clarity. This griding 
approach provides the best compromise. A schematic of a cross section of the internal 
nozzle is shown in Fig. 6.6.36. The band of polar grid properly meets the external grid 
with C 1 continuity at the exit of the afterbody, as well as provides the normal distance 
from the solid interior walls for the length scale necessary for the algebraic turbulence 
model. Plus, having the polar grid interface with the H-H grid reduces the difficulty in 
maintaining the polar lines normal to the solid surfaces. The polar grid meshed with 
the H-H grid quite easily, with the caveat that the cell sizes across the interface had to 
be within 20% of each other. This is the same rate of change that is allowed between 
cells in the absence of an interface. If this were not enforced there would be metric 
discontinuities across the interface, which would contaminate the solution. 
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Figure 6.6.33 Afterbody Surface and Exterior Polar Grid Configuration. 
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Figure 6.6.34 Afterbody Grid Geometry. 
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Figure 6.6.35 Internal Nozzle with Combined H-H and Polar Grid Topologies. 
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Figure 6.6.36 Schematic of H-H Grid and Polar Grid Topologies Interfacing. 

The first case tested was for M 0 Q - 0.60, a = 0.0, M jet = 0.33, p je t/Poc = 3 7 1, 
Tjet/Too = 0.97, and R = 273,000/(unit length), with the afterbody being 63.04 units 
long. The grid used in Ref [105] was obtained, which contained only the external 
geometry. This grid was used as a preliminary grid only, because it is not well suited 
for multigrid applications since locations of geometrical change were not at multigridable 
indices. This causes the physical locations to change for different multigrid levels, which 
in turn can introduce error. Calculations were performed on this grid, with the assumption 
that the plume emanating from the nozzle would remain the same size as the cross 
sectional area of the nozzle. This grid configuration assumed the end of the afterbody 
had a sharp trailing edge, rather than incorporate the flat surface which separates the 
external surface from the internal nozzle. The difference is shown schematically in 
Fig. 6.6.37. The grid had 129x33x65 points, with 129 points streamwise, 33 points 
distributed circumferentially, and 65 points normal to the exterior surface. The cell 
spacing normal to the surface was lxlO -4 units. Only Mlf of the body was grided. 
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because the angle of attack and yaw were zero; therefore the flow was symmetric about 
the longitudinal axis. The results for this case are shown in Fig. 6.6.38, for the C v data 
on the top of the boat-tail, and Fig. 6.6.39 for the C p data on the side of the boat-tail. 
Good agreement was obtained for both locations in comparison with the experimental 
data [105], with the exception of the end of the boat-tail. These results are comparable to 
the numerical results of Compton, Thomas, Abeyounis, and Mason [105], for the same 
case where they used a Baldwin-Lomax algebraic turbulence model. The same grid was 
used for a case at M ^ = 0.80, a = 0.0, Mj e t = 0.33, Pjet/Po o = 3.71, Tjet/Too — 0.99, 
and R = 309, 000/ (unit length). The results for the top and side wall C v data are shown 
in Figs. 6.6.40 and 6.6.41, respectively. Again, the agreement with the experimental 
data, was good [105], and comparable to the numerical results of Compton et al [105]. 


Preliminary Configuration 
Geometry 



* 



Figure 6.6.37 Comparison of Preliminary Configuration 
Geometry and Actual Configuration Geometry. 
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Figure 6.6.38 Preliminary Afterbody Nozzle Top Wall Pressure Coefficient. 



Figure 6.6.39 Preliminary Afterbody Nozzle Side Wall Pressure Coefficient. 






Figure 6.6.40 Preliminary Afterbody Nozzle Top Wall Pressure Coefficient 



Figure 6.6.41 Preliminary Afterbody Nozzle Top Wall Pressure Coefficient 





Since the interior grid needed to be computed, along with the concerns about the 
external grid being ill suited for multigrid, a new multigrid compatible grid was generated. 
This grid had 177 points in the streamwise direction, with 97 points on the body. Since 
the flow was to be symmetric about the longitudinal axis, only the quarter plane was 
generated, with 33 points in the circumferential direction, and 49 points in the direction 
normal to the body. The cell spacing normal to the body on the boat-tail was 5x10 units. 

The first task after generating the new grid was to make a comparison between it and 
the preliminary grid, using only the external grid configurations. The results for the top 
and side wall C p data for the case of = 0.60, a = 0.0, M jet = 0.33, p je t/p<x = 3.71, 
Tjet/Too = 0-97, and R = 273,000/(unit length), are shown in Figs. 6.6.42 and 6.6.43, 
respectively. As can be seen, the comparison between the two different grids is good. 
The new grid was then tested for the case of M <*> = 0.80, a = 0.0, M jet = 0.33, 
Pjet / Po o = 3.71, Tjet/Too = 0.99, and R = 309,000/(unit length). Again the agreement 
was good, as can be seen by comparisons of the C v data for the top and side walls in 
Figures 6.6.44 and 6.6.45, respectively. For these cases the Baldwin-Lomax algebraic 
turbulence model was active only over the solid wall surfaces, and the shear line, which 
should have been between the plume and the free-steam flow was treated as an inviscid 

surface. 
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Figure 6.6.45 Comparison of Single-Block 
Afterbody Grids for Side Wall Pressure Coefficient 
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After verifying that the newly generated grid performed adequately, the test case 
of Moo = 0.60, q = 0.0, M )ei = 0.33, p je tlPoo = 3.71, T^t/Too = 0.97, and 
R — 273, 000/ (unit length) was tested again, utilizing the grid that included the thickness 
between the exterior surface and the internal grid at the end of the afterbody, as shown 
in Fig. 6.6.46. As can be seen by the comparison of C p predictions for the top and 
side walls in Fig. 6.6.47 and 6.6.48, the agreement with the experimental data [105] was 
very similar to that of the single-block configuration. However, the blunt geometry was 
somewhat better in predicting the aft portion of the sidewall pressure coefficient. 



Figure 6.6.46 Schematic of Afterbody for Two Block External Configuration. 
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Figure 6.6.48 Comparison of Single-Block and Two-Block 
Afterbody Nozzle Side Wall Pressure Coefficient. 
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The same case was examined again employing both the plume grids and the internal 
nozzle grids, as shown in Fig. 6.6.49. This gave a total of four blocks. The external 
block extended only to the end of the afterbody, as did the interior polar block. The H-H 
grid block began in the circular region of the internal nozzle, as did the interior polar 
grid, and extended all the way to the exit plane of the numerical domain, which was 
approximately 1 10 units downstream from the end of the afterbody. The fourth block 



Figure 6.6.49 Schematic of Four-Block Internal and External Afterbody Configuration 

began at the end of the afterbody, where that face of the block incorporated three boundary 
conditions. One was to interface with the polar grid block from the interior nozzle. The 
second boundary condition was to represent the solid wall thickness between the internal 
and external surfaces of the afterbody, and the final condition was to interface with the 
external block of the afterbody. Again, the algebraic turbulence model was employed 
only on the solid wall surfaces of the afterbody, including both the interior and exterior 
surfaces. The results of this case are shown as comparisons of C p data for the top and 


137 


side walls in Fig. 6.6.50 and 6.6.51, where there is good agreement with the experimental 
data [105] and the numerical results of the two-block case, except at the trailing edge of 
the boat-tail. One reason the four-block configuration predicted higher pressures on the 
surface of the boat-tail maybe because as the flow exits the nozzle, it has a velocity that 
is not parallel to the horizontal axis, as was the assumption for the one- and two-block 
configurations. There isaw velocity component due to the divergence of the nozzle, 
which forces the shear layer to turn upward, and generates a higher pressure upstream, 
on the boat-tail. Also, the flow coming out of the nozzle is supersonic, so it is going 
to expand if the surrounding pressure regions allow it. This expansion will also cause 
the shear line to turn outward. Flow from on the boat-tail definitely expands into the 
thickness region at the end of the afterbody, and when subsonic flow expands its velocity 
drops and its pressure increases. This higher pressure in the thickness region can influence 
the pressures upstream on the boat-tail, and there is no direct pressure correction term 
in the Baldwin-Lomax algebraic turbulence model; therefore the higher pressures at the 
end of the afterbody, in the thickness region, may erroneously influence the pressures on 
the boat-tail. Also in this region, there appears to be a re-circulation bubble which can 
add to the higher upstream pressures. 
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Figure 6.6.50 Comparison of Two-Block and Four-Block 
Afterbody Nozzle Top Wall Pressure Coefficient. 



Figure 6.6.5 1 Comparison of Two-Block and Four-Block 
Afterbody Nozzle Side Wall Pressure Coefficient. 
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The same four-block configuration was tested for the Moo - 0.80, a • , 
Mj„ = 0.33, Pl e,IPoo = 3.71, T, ,,/Too = 0.99, and R = 309, 000/ (unit length) 
case. Again, the pressures on the boat-tail were higher than what the single-block 
configuration predicted, as is displayed by the C r data for the top and side walls in 

Figs. 6.6.52 and 6.6.53. 



i/c 


Figure 6.6.52 Comparison of Single-Block and Four-Block 
Afterbody Nozzle Top Wall Pressure Coefficient. 
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Figure 6.6.53 Comparison of Single-Block and Four-Block 
Afterbody Nozzle Side Wall Pressure Coefficient. 
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Chapter 7 CONCLUSIONS 


This work was aimed at developing a state of the art computer code, capable of 
evaluating arbitrary, complex geometric configurations. The computer code developed 
utilizes upwind solvers, van Leer’s flux-vector splitting and Roe s flux-difference splitting, 
coupled with an explicit modified Runge-Kutta time integration scheme. To accelerate 
the rate of convergence the techniques of FAS multigrid, variable coefficient residual 
smoothing and constant coefficient corrector smoothing were incorporated. Development 
of the computer code was an attempt to capitalize on the benefits of both the upwind 
implicit time integration and the central-difference explicit Runge-Kutta time integration 
methods. Although the upwind methods require more CPU time for the flux evaluations, 
they provide better shock capturing and less’ smearing of contact discontinuities than 
the central-difference methods. However, the multistage central difference methods 
are generally more capable of smoothing out the high frequency errors, allowing the 
multigrid acceleration techniques to restrict information without causing aliasing. For 
the upwind computer codes to provide a sufficient amount of smoothing an implicit 
time integration is usually used, which again requires more CPU time. The present 
computer code incorporated a modified Runge-Kutta time integration technique, very 
similar to what the central-difference codes employ, in an effort to attain the speed 
of the central-difference codes, coupled with the accuracy of the upwind codes. The 
thrust of this work was to aid in the development of the computer code and provide 
the flexibility of multi-block capabilities. It was desired that this computer code be 
capable of accommodating any configuration that could be defined by six-sided blocks. To 
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incorporate multigrid techniques efficiently requires an understanding of how computers 
store data, and operate. The same can be said for incorporating multi-block capabilities. 
Many of the qualities a computer program requires for multigrid techniques are the same 
for multi-block capabilities, which are mainly grid level, block, and topology independent 
subroutines. These attributes allow for having only one subroutine to execute a particular 
task, regardless of what grid level, block, or mesh topology is being evaluated. This 
leads to having a generic computer program, which allows it will be applicable to any 
configuration. 

The test cases were chosen based on the criteria of validating the computer code and 
displaying its versatility. The corner flow required three-dimensional flow capabilities. 
Also, when it was divided into eight blocks, it demanded special interface requirements 
allowing the multigrid prolongation calculations to be independent of the interface 
locations. This type of interface condition would be common at the back end of an 
aircraft engine, where internal and external flows meet, or where flows join from around 
different solid wall body parts. Good agreement with other numerical results [98, 99] 
were obtained. This case was performed using multigrid acceleration, with and without 
variable coefficient residual smoothing. Since this case was supersonic flow, and the grids 
were not highly stretched, the residual smoothing did not aid the rate of convergence. 
This acceleration technique is best suited for stretched grids, such as viscous grids, and 
obviously more compatible with subsonic and transonic flows, than for supersonic flows. 

The pseudo two-dimensional jet plume required multiple inlet boundary conditions, 
and provided a testing ground for comparing how well different solvers would compute 
the slip line and shock surface. Favorable comparisons were made with the validated 
shock fitting code of Salas [101]. For the upwind methods. Roe’s flux-difference splitting 
and van Leer’s flux-vector splitting, their ability to accurately predict shocks and slip lines 


143 



was dependent on how these structures were aligned with the computational grid. The 
flux-difference splitting method should provide better slip line resolution than flux-vector 
splitting, because it is an approximate Riemann solver, which by design is capable of 
resolving contact surfaces. The flux-vector splitting method resolved the slip line just 
as accurately as the flux-difference splitting method, partly because the slip line was not 
directly aligned with the grid. If it was better aligned, the flux-difference splitting method 
should have provided better resolution. Based on this case, it was determined that both 
upwind methods would provide acceptable results for the rest of the cases studied. 

The laminar and turbulent flat plate cases were computed with great success. The 
laminar results compared well with the Blasius solution, indicating that the viscous terms 
had been properly incorporated into the governing equations. The turbulent results 
compared well with the analytical velocity profile and skin friction, plus the law of 
the wall plot. Thus indicating that the Baldwin-Lomax algebraic turbulence model [72] 
was functioning properly for attached turbulent flows. Full multigrid acceleration was 
used, with the cell weighted residual smoothing. These two techniques made a large 
improvement in the convergence rate for these two cases. 

The ONERA M6 wing cases were true three-dimensional attached turbulent flows. 
Results for both cases compared well with the experimental data [102] and other numerical 
results [86, 103]. These cases were computed using FMG acceleration and variable 
coefficient residual smoothing. Again, these acceleration techniques provided improved 
convergence rates. This configuration also placed special requirements on the interface 
conditions, which are needed when a block face has an interface and a wall boundary 
condition next to each other. 

The final case studied was the afterbody configuration. This configuration required 
a multi-block code in order to accommodate the internal and external geometries. The 
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results for this case are preliminary due to the quality of the grid. The resolution on the 
solid surfaces was sufficient for the formulation of the boundary layers, but at locations 
of large geometrical change many cells became abnormally skewed. Much could be done 
to improve the grid, with the proper experience and software. The afterbody grid went 
through four revisions before a final selection was made. Three main areas required 
improvements: the nose of the afterbody, the sharp compression in the nozzle geometry 
(which started just upstream of the throat) and the thickness region at the end of the 
afterbody. The nose created numerical difficulty because of the polar axis coupled with 
viscous compatible cell spacing. Relaxing the cell spacing normal to the nose and along 
the polar axis greatly relieved some of the numerical difficulty, allowing this region 
to numerically converge. The thickness region, at the end of the afterbody, required 
a smooth continuous variation of cell spacing from the interior nozzle to the exterior 
surface. This was a difficult procedure, especially at the corners, because the grid lines 
from the interior nozzle had to match the grid lines of the external surface. Also, the 
thickness of this region, at the end of the afterbody, changed at the comers of the boat- 
tail. The top and bottom thicknesses were 0.0395 units, and the sides were 0.125 units. 
This made having smooth variations of cell sizes from the interior nozzle to the exterior 
surface more difficult. The start of the compression region for the interior nozzle requires 
true three-dimensional grid surfaces to reduce the cell skewness. This region is where 
many of the numerical difficulties remained. Much more sophisticated grid generation 
techniques would be required to effectively approach this problem. 

This case required a multi-block computer program to simultaneously accommodate 
the internal and external domains. The computer program performed well based on the 
pressure coefficient results obtained for the two flow cases. Further improvements in the 
grid quality may not effect these results, but employing a different type of turbulence 
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model definitely could, especially if it retains the turbulent kinetic energy term, which 
is eliminated in algebraic turbulence models. Note that if the turbulent kinetic energy 
term is actively accounted for, then there will be a direct link with the static pressure 
through the equation of state [106]. 

Overall, the computer code performed well. Its design enabled it to handle many dif- 
ferent configurations without requiring any alteration of the source code. This flexibility 
makes it a very versatile tool in examining many types of flow configurations. Also, it 
provided competitive solutions when compared with other numerical results [86, 98, 99, 
101, 103] and experimental data [102]. The multigrid acceleration decreased the amount 
of work units required to obtain a solution for all cases, except supersonic Euler flow on 
an unstretched grid. Since the multigrid process executes averaging that resembles ellip- 
tic information propagation, it sent information upstream in the supersonic flow cases, 
which can obviously reduce the convergence rate. A typical example of this was the 
supersonic corner flow. As for subsonic and transonic flows, and flows on viscous grids, 
the multigrid acceleration was a great aid to accelerate the rate of convergence. The vari- 
able coefficient residual smoothing had the best impact on the viscous grids as well. This 
is due primarily to the fact that the residual smoothing coefficients are generated relative 
to a specified value of cell aspect ratios and spectral radii. This smoothing technique did 
aid in the rate of convergence for viscous flow cases. The corrector smoothing, which 
used constant coefficients, never accelerated the rate of convergence for the cases tested. 
It also was intended to aid the viscous flow cases, but for the range of Mach numbers and 
grids employed in the work presented, it actually reduced the rate of convergence for test 
performances, and therefore was not used in the final analysis of the flow configurations. 

If the numerical efficiency of the computer program was measured based on the 
amount of work units necessary to obtain global flow values, the computer code did 
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well. Generally, it required approximately 150 to 200 work units of FMG, using the pure 
second order upwind extrapolation, four stage modified Runge-Kutta time integration, 
coupled with variable coefficient residual smoothing, for the afterbody cases. The global 
residual, which was based on the change of density, did not converge very well. Is was 
case dependent as well as being very grid dependent, which is why four different grids 
were generated for the afterbody. The wing cases converged about three orders under 
the same execution conditions. Again, the global quantities were obtained with just a 
couple hundred work units of FMG. The flat plate cases converged better, based on the 
amount the residual dropped, but it did require more work units until the entire boundary 
layer was completely developed, because of the high number of cells in the boundary 
layer. All of the inviscid cases converged well. Unfortunately the more numerically 
demanding cases, such as the viscous cases, require more frequency damping than is 
currently being provided by the modified Runge-Kutta time integration scheme. Much 
more attention needs to be directed toward modeling the frequencies that are being 
produced with this type of numerical approach, so that a more accurate set of modified 
Runge-Kutta coefficients can be determined. This problem has been investigated by a 
number of scientist with some success, but much still needs to be done. One issue 
is the incorporation of such acceleration techniques as multigrid, residual and corrector 
smoothing. An other issue is to actually incorporate the van Leer’s flux-vector splitting 
and Roe’s flux-difference splitting techniques in the model equations. 

True multigrid performance is rarely seen in any of the complicated flow configura- 
tions, for any method, but especially for upwind methods. Having a computer program 
that is completely independent of cell aspect ratios is still a goal that has not been 
achieved by anyone in the CFD community. 
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Appendix A Full Navier-Stokes Equations 
in Body-Fitted Coordinates 

For the ^-direction 

F v { 1) = [0.0] (A.l) 


«<(&* + tl + t 2 z) + Unihit* + Vyty + + 

A* rU w <(|Cx^r + + Cz6) + 
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For the //-direction 
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For the (-direction 
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Appendix B Multigrid Restriction and 
Prolongation Operations 

B.l Restriction Operation 

The best way to explain the restriction procedure is by showing the operation for a 
two-dimensional case. Starting with the fine grid, which is comprised of the small cells 
whose average values are designated by the empty circles, the coarse grid is generated 
by removing every other fine grid line (see Fig. B.l). Thus, producing the larger cells 
which are designated by the thicker borders. The shaded octagons, centered in the larger 
cells, represent the coarser grid cells’ centered values. As shown in Fig. B.l, a volume 
weighted averaging of the flow values from the fine grid cells is used to provide coarse 
grid values, which are to represent the solution of the fine grid on the coarser grid. 



Figure B.l Two-Dimensional Restriction Operation. 


161 




In three-dimensional space, the same type of averaging of the fine grid values is 
performed on the two fine grid planes that surround each coarse grid plane. To show 
this, it is necessary to indicate where the coarse grid and fine grid cell centers are located 
geometrically. This can be seen in Fig. B.2, which represents a three-dimensional 
volume. The fine grid is indicated by the thin lines, and the coarse grid by the thicker 
lines. The fine grid cell centers are at the centers of all the small cells, and the coarse grid 
cell centers are at the centers of the larger cells, indicated by the thicker lines. Cutting 
some of the cells away and putting in four planes, two representing fine grid cell centers 
and two representing coarse grid cell centers, shows the spacial relation of the coarse 
and fine grid cell centers, as can be seen in Fig. B.1.3. This figure identifies the various 
planes that are needed for both restriction and prolongation operations. Each coarse grid 
cell centered plane has a fine grid cell centered plane on both sides of it. Focusing just 
on these three planes, as shown in Fig. B.1.4, this averaging will involve eight fine 
grid cells to produce one coarse grid cell value. The volume weighted averaging uses 
the actual physical cell’s volumes, and not those of the computational domain, which 
are unit volumes. 
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Figure B.2 Three-Dimensional Fine and Coarse Grid. 
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Figure B.1.3 Three-Dimensional Fine and Coarse Grid Cell Centers 
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Figure B.1.4 Three-Dimensional Restriction Operation 


B.2 Prolongation Operation 

The prolongation operation is done in computational space. For the two-dimensional 
case a bilinear interpolation is performed among the coarse grid cells, which are indicated 
by the shaded octagons, shown in Fig. B.2.5. Each set of four shaded octagons has four 
empty circles within the octagons’ perimeter. These circles represent the fine grid cells. 
The octagon that is closest to the circle has the largest weight factor, which gives it the 
most influence on the circle’s value. How the prolongation values, for each circle in the 
perimeter, are obtained is indicated in the four different patterns shown in Fig. B.2.5. 
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Figure B.2.5 Two-Dimensional Prolongation Operation. 


In three-dimensional space, the prolongation operation is again done in the computa- 
tional domain using a tri-linear interpolation. Referring back to Fig. B.1.3, the five planes 
are re-drawn in Fig. B.2.6, with circles added to indicate the cell centers. First a bilinear 
interpolation is performed on each coarse grid cell center plane, as shown in Fig. B.2.6. 
Then a linear interpolation between the two coarse grid prolongations is performed to 
give the value for the specified fine grid cell. As can be seen in Fig. B.2.6, there is a set 
of eight coarse grid cells that provide information to produce eight fine grid cells. 
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Figure B.2.6 Three-Dimensional Prolongation Operation. 


167 










